<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<title>MC-GPU: MC-GPU_kernel_v1.3.cu Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<script type="text/javascript">
  $(document).ready(initResizable);
</script>


</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">MC-GPU
   
   </div>
   
  </td>
  
  
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.6.1 -->
  <div id="navrow1" class="tabs">
    <ul class="tablist">
      <li><a href="index.html"><span>Main&#160;Page</span></a></li>
      <li><a href="annotated.html"><span>Classes</span></a></li>
      <li class="current"><a href="files.html"><span>Files</span></a></li>
    </ul>
  </div>
  <div id="navrow2" class="tabs2">
    <ul class="tablist">
      <li><a href="files.html"><span>File&#160;List</span></a></li>
      <li><a href="globals.html"><span>File&#160;Members</span></a></li>
    </ul>
  </div>
</div>
<div id="side-nav" class="ui-resizable side-nav-resizable">
  <div id="nav-tree">
    <div id="nav-tree-contents">
    </div>
  </div>
  <div id="splitbar" style="-moz-user-select:none;" 
       class="ui-resizable-handle">
  </div>
</div>
<script type="text/javascript">
  initNavTree('MC-GPU__kernel__v1_83_8cu.html','');
</script>
<div id="doc-content">
<div class="header">
  <div class="headertitle">
<div class="title">MC-GPU_kernel_v1.3.cu</div>  </div>
</div><!--header-->
<div class="contents">
<a href="MC-GPU__kernel__v1_83_8cu.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment"></span>
<a name="l00002"></a>00002 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00003"></a>00003 <span class="comment"></span><span class="comment">//</span>
<a name="l00004"></a>00004 <span class="comment">//              ****************************</span>
<a name="l00005"></a>00005 <span class="comment">//              *** MC-GPU , version 1.3 ***</span>
<a name="l00006"></a>00006 <span class="comment">//              ****************************</span>
<a name="l00007"></a>00007 <span class="comment">//                                          </span><span class="comment"></span>
<a name="l00008"></a>00008 <span class="comment">//!  Definition of the CUDA GPU kernel for the simulation of x ray tracks in a voxelized geometry.</span>
<a name="l00009"></a>00009 <span class="comment">//!  This kernel has been optimized to yield a good performance in the GPU but can still be</span>
<a name="l00010"></a>00010 <span class="comment">//!  compiled in the CPU without problems. All the CUDA especific commands are enclosed in</span>
<a name="l00011"></a>00011 <span class="comment">//!  pre-processor directives that are skipped if the parameter &quot;USING_CUDA&quot; is not defined</span>
<a name="l00012"></a>00012 <span class="comment">//!  at compilation time.</span>
<a name="l00013"></a>00013 <span class="comment"></span><span class="comment">//</span>
<a name="l00014"></a>00014 <span class="comment">//        ** DISCLAIMER **</span>
<a name="l00015"></a>00015 <span class="comment">//</span>
<a name="l00016"></a>00016 <span class="comment">// This software and documentation (the &quot;Software&quot;) were developed at the Food and</span>
<a name="l00017"></a>00017 <span class="comment">// Drug Administration (FDA) by employees of the Federal Government in the course</span>
<a name="l00018"></a>00018 <span class="comment">// of their official duties. Pursuant to Title 17, Section 105 of the United States</span>
<a name="l00019"></a>00019 <span class="comment">// Code, this work is not subject to copyright protection and is in the public</span>
<a name="l00020"></a>00020 <span class="comment">// domain. Permission is hereby granted, free of charge, to any person obtaining a</span>
<a name="l00021"></a>00021 <span class="comment">// copy of the Software, to deal in the Software without restriction, including</span>
<a name="l00022"></a>00022 <span class="comment">// without limitation the rights to use, copy, modify, merge, publish, distribute,</span>
<a name="l00023"></a>00023 <span class="comment">// sublicense, or sell copies of the Software or derivatives, and to permit persons</span>
<a name="l00024"></a>00024 <span class="comment">// to whom the Software is furnished to do so. FDA assumes no responsibility</span>
<a name="l00025"></a>00025 <span class="comment">// whatsoever for use by other parties of the Software, its source code,</span>
<a name="l00026"></a>00026 <span class="comment">// documentation or compiled executables, and makes no guarantees, expressed or</span>
<a name="l00027"></a>00027 <span class="comment">// implied, about its quality, reliability, or any other characteristic. Further,</span>
<a name="l00028"></a>00028 <span class="comment">// use of this code in no way implies endorsement by the FDA or confers any</span>
<a name="l00029"></a>00029 <span class="comment">// advantage in regulatory decisions.  Although this software can be redistributed</span>
<a name="l00030"></a>00030 <span class="comment">// and/or modified freely, we ask that any derivative works bear some notice that</span>
<a name="l00031"></a>00031 <span class="comment">// they are derived from it, and any modified versions bear some notice that they</span>
<a name="l00032"></a>00032 <span class="comment">// have been modified.</span>
<a name="l00033"></a>00033 <span class="comment">//                                                                            </span>
<a name="l00034"></a>00034 <span class="comment">//</span><span class="comment"></span>
<a name="l00035"></a>00035 <span class="comment">//!                     @file    MC-GPU_kernel_v1.3.cu</span>
<a name="l00036"></a>00036 <span class="comment">//!                     @author  Andreu Badal (Andreu.Badal-Soler@fda.hhs.gov)</span>
<a name="l00037"></a>00037 <span class="comment">//!                     @date    2012/12/12</span>
<a name="l00038"></a>00038 <span class="comment"></span><span class="comment">//                       -- Original code started on:  2009/04/14</span>
<a name="l00039"></a>00039 <span class="comment">//</span><span class="comment"></span>
<a name="l00040"></a>00040 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00041"></a>00041 <span class="comment"></span>
<a name="l00042"></a>00042 
<a name="l00043"></a>00043 <span class="comment"></span>
<a name="l00044"></a>00044 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00045"></a>00045 <span class="comment">//!  Initialize the image array, ie, set all pixels to zero</span>
<a name="l00046"></a>00046 <span class="comment">//!  Essentially, this function has the same effect as the command: </span>
<a name="l00047"></a>00047 <span class="comment">//!   &quot;cutilSafeCall(cudaMemcpy(image_device, image, image_bytes, cudaMemcpyHostToDevice))&quot;;</span>
<a name="l00048"></a>00048 <span class="comment">//!  </span>
<a name="l00049"></a>00049 <span class="comment">//!  CUDA performs some initialization work the first time a GPU kernel is called.</span>
<a name="l00050"></a>00050 <span class="comment">//!  Therefore, calling a short kernel before the real particle tracking is performed</span>
<a name="l00051"></a>00051 <span class="comment">//!  may improve the accuracy of the timing measurements in the relevant kernel.</span>
<a name="l00052"></a>00052 <span class="comment">//!  </span>
<a name="l00053"></a>00053 <span class="comment">//!       @param[in,out] image   Pointer to the image array.</span>
<a name="l00054"></a>00054 <span class="comment">//!       @param[in] pixels_per_image  Number of pixels in the image (ie, elements in the array).</span>
<a name="l00055"></a>00055 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00056"></a>00056 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00057"></a>00057 <span class="preprocessor"></span>__global__
<a name="l00058"></a>00058 <span class="keywordtype">void</span> init_image_array_GPU(<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span>* image, <span class="keywordtype">int</span> pixels_per_image)
<a name="l00059"></a>00059 {
<a name="l00060"></a>00060   <span class="keywordtype">int</span> my_pixel = threadIdx.x + blockIdx.x*blockDim.x;
<a name="l00061"></a>00061   <span class="keywordflow">if</span> (my_pixel &lt; pixels_per_image)
<a name="l00062"></a>00062   {
<a name="l00063"></a>00063     <span class="comment">// -- Set the current pixel to 0 and return, avoiding overflow when more threads than pixels are used:</span>
<a name="l00064"></a>00064     image[my_pixel] = (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)(0);    <span class="comment">// Initialize non-scatter image</span>
<a name="l00065"></a>00065     my_pixel += pixels_per_image;                     <span class="comment">//  (advance to next image)</span>
<a name="l00066"></a>00066     image[my_pixel] = (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)(0);    <span class="comment">// Initialize Compton image</span>
<a name="l00067"></a>00067     my_pixel += pixels_per_image;                     <span class="comment">//  (advance to next image)</span>
<a name="l00068"></a>00068     image[my_pixel] = (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)(0);    <span class="comment">// Initialize Rayleigh image</span>
<a name="l00069"></a>00069     my_pixel += pixels_per_image;                     <span class="comment">//  (advance to next image)</span>
<a name="l00070"></a>00070     image[my_pixel] = (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)(0);    <span class="comment">// Initialize multi-scatter image</span>
<a name="l00071"></a>00071   }
<a name="l00072"></a>00072 }
<a name="l00073"></a>00073 
<a name="l00074"></a>00074 <span class="comment">// ////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00075"></a>00075 <span class="comment">// //!  Initialize the dose deposition array, ie, set all voxel doses to zero</span>
<a name="l00076"></a>00076 <span class="comment">// //!  </span>
<a name="l00077"></a>00077 <span class="comment">// //!       @param[in,out] dose   Pointer to the dose mean and sigma arrays.</span>
<a name="l00078"></a>00078 <span class="comment">// //!       @param[in] num_voxels_dose  Number of voxels in the dose ROI (ie, elements in the arrays).</span>
<a name="l00079"></a>00079 <span class="comment">// ////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00080"></a>00080 <span class="comment">// __global__</span>
<a name="l00081"></a>00081 <span class="comment">// void init_dose_array_GPU(ulonglong2* voxels_Edep, int num_voxels_dose)</span>
<a name="l00082"></a>00082 <span class="comment">// {  </span>
<a name="l00083"></a>00083 <span class="comment">//   int my_voxel = threadIdx.x + blockIdx.x*blockDim.x;</span>
<a name="l00084"></a>00084 <span class="comment">//   register ulonglong2 ulonglong2_zero;</span>
<a name="l00085"></a>00085 <span class="comment">//   ulonglong2_zero.x = ulonglong2_zero.y = (unsigned long long int) 0;</span>
<a name="l00086"></a>00086 <span class="comment">//   if (my_voxel &lt; num_voxels_dose)</span>
<a name="l00087"></a>00087 <span class="comment">//   {</span>
<a name="l00088"></a>00088 <span class="comment">//     dose[my_voxel] = ulonglong2_zero;    // Set the current voxel to (0,0) and return, avoiding overflow</span>
<a name="l00089"></a>00089 <span class="comment">//   }</span>
<a name="l00090"></a>00090 <span class="comment">// }</span>
<a name="l00091"></a>00091 
<a name="l00092"></a>00092 <span class="preprocessor">#endif</span>
<a name="l00093"></a>00093 <span class="preprocessor"></span>
<a name="l00094"></a>00094  <span class="comment"></span>
<a name="l00095"></a>00095 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00096"></a>00096 <span class="comment">//!  Main function to simulate x-ray tracks inside a voxelized geometry.</span>
<a name="l00097"></a>00097 <span class="comment">//!  Secondary electrons are not simulated (in photoelectric and Compton </span>
<a name="l00098"></a>00098 <span class="comment">//!  events the energy is locally deposited).</span>
<a name="l00099"></a>00099 <span class="comment">//!</span>
<a name="l00100"></a>00100 <span class="comment">//!  The following global variables, in  the GPU __constant__ memory are used:</span>
<a name="l00101"></a>00101 <span class="comment">//!           voxel_data_CONST, </span>
<a name="l00102"></a>00102 <span class="comment">//!           source_energy_data_CONST,</span>
<a name="l00103"></a>00103 <span class="comment">//!           detector_data_CONST, </span>
<a name="l00104"></a>00104 <span class="comment">//!           mfp_table_data_CONST.</span>
<a name="l00105"></a>00105 <span class="comment">//!</span>
<a name="l00106"></a>00106 <span class="comment">//!       @param[in] history_batch  Particle batch number (only used in the CPU version when CUDA is disabled!, the GPU uses the built-in variable threadIdx)</span>
<a name="l00107"></a>00107 <span class="comment">//!       @param[in] num_p  Projection number in the CT simulation. This variable defines a specific angle and the corresponding source and detector will be used.</span>
<a name="l00108"></a>00108 <span class="comment">//!       @param[in] histories_per_thread   Number of histories to simulate for each call to this function (ie, for GPU thread).</span>
<a name="l00109"></a>00109 <span class="comment">//!       @param[in] seed_input   Random number generator seed (the same seed is used to initialize the two MLCGs of RANECU).</span>
<a name="l00110"></a>00110 <span class="comment">//!       @param[in] voxel_mat_dens   Pointer to the voxel densities and material vector (the voxelized geometry), stored in GPU glbal memory.</span>
<a name="l00111"></a>00111 <span class="comment">//!       @param[in] mfp_Woodcock_table    Two parameter table for the linear interpolation of the Woodcock mean free path (MFP) (stored in GPU global memory).</span>
<a name="l00112"></a>00112 <span class="comment">//!       @param[in] mfp_table_a   First element for the linear interpolation of the interaction mean free paths (stored in GPU global memory).</span>
<a name="l00113"></a>00113 <span class="comment">//!       @param[in] mfp_table_b   Second element for the linear interpolation of the interaction mean free paths (stored in GPU global memory).</span>
<a name="l00114"></a>00114 <span class="comment">//!       @param[in] rayleigh_table   Pointer to the table with the data required by the Rayleigh interaction sampling, stored in GPU global memory.</span>
<a name="l00115"></a>00115 <span class="comment">//!       @param[in] compton_table   Pointer to the table with the data required by the Compton interaction sampling, stored in GPU global memory.</span>
<a name="l00116"></a>00116 <span class="comment">//!       @param[in,out] image   Pointer to the image vector in the GPU global memory.</span>
<a name="l00117"></a>00117 <span class="comment">//!       @param[in,out] dose   Pointer to the array containing the 3D voxel dose (and its uncertainty) in the GPU global memory.</span>
<a name="l00118"></a>00118 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00119"></a>00119 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00120"></a>00120 <span class="preprocessor"></span>__global__ <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ab1e2e311c809fe74d093498e53c9c954" title="Initialize the image array, ie, set all pixels to zero Essentially, this function has the same effect...">track_particles</a>(<span class="keywordtype">int</span> histories_per_thread,
<a name="l00121"></a>00121                                 <span class="keywordtype">int</span> num_p,      <span class="comment">// For a CT simulation: allocate space for up to MAX_NUM_PROJECTIONS projections.</span>
<a name="l00122"></a>00122                                 <span class="keywordtype">int</span> seed_input,
<a name="l00123"></a>00123                                 <span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span>* image,
<a name="l00124"></a>00124                                 <a class="code" href="structulonglong2.html">ulonglong2</a>* voxels_Edep,
<a name="l00125"></a>00125                                 <a class="code" href="structfloat2.html">float2</a>* voxel_mat_dens,
<a name="l00126"></a>00126                                 <a class="code" href="structfloat2.html">float2</a>* mfp_Woodcock_table,
<a name="l00127"></a>00127                                 <a class="code" href="structfloat3.html">float3</a>* mfp_table_a,
<a name="l00128"></a>00128                                 <a class="code" href="structfloat3.html">float3</a>* mfp_table_b,
<a name="l00129"></a>00129                                 <span class="keyword">struct</span> <a class="code" href="structrayleigh__struct.html" title="Structure storing the data of the Rayleigh interaction sampling model (equivalent to PENELOPE&#39;s commo...">rayleigh_struct</a>* rayleigh_table,
<a name="l00130"></a>00130                                 <span class="keyword">struct</span> <a class="code" href="structcompton__struct.html" title="Structure storing the data of the Compton interaction sampling model (equivalent to PENELOPE&#39;s common...">compton_struct</a>* compton_table,
<a name="l00131"></a>00131                                 <span class="keyword">struct</span> <a class="code" href="structdetector__struct.html" title="Structure storing the data defining the x-ray detector.">detector_struct</a>* detector_data_array,
<a name="l00132"></a>00132                                 <span class="keyword">struct</span> <a class="code" href="structsource__struct.html" title="Structure storing the data defining the source model (except for the energy spectrum).">source_struct</a>* source_data_array, 
<a name="l00133"></a>00133                                 <a class="code" href="structulonglong2.html">ulonglong2</a>* materials_dose)
<a name="l00134"></a>00134 <span class="preprocessor">#else</span>
<a name="l00135"></a><a class="code" href="MC-GPU__v1_83_8h.html#ab1e2e311c809fe74d093498e53c9c954">00135</a> <span class="preprocessor"></span>           <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ab1e2e311c809fe74d093498e53c9c954" title="Initialize the image array, ie, set all pixels to zero Essentially, this function has the same effect...">track_particles</a>(<span class="keywordtype">int</span> history_batch,             <span class="comment">// This variable is not required in the GPU, it uses the thread ID           </span>
<a name="l00136"></a>00136                                 <span class="keywordtype">int</span> histories_per_thread,
<a name="l00137"></a>00137                                 <span class="keywordtype">int</span> num_p,
<a name="l00138"></a>00138                                 <span class="keywordtype">int</span> seed_input,
<a name="l00139"></a>00139                                 <span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span>* image,
<a name="l00140"></a>00140                                 <a class="code" href="structulonglong2.html">ulonglong2</a>* voxels_Edep,
<a name="l00141"></a>00141                                 <a class="code" href="structfloat2.html">float2</a>* voxel_mat_dens,
<a name="l00142"></a>00142                                 <a class="code" href="structfloat2.html">float2</a>* mfp_Woodcock_table,
<a name="l00143"></a>00143                                 <a class="code" href="structfloat3.html">float3</a>* mfp_table_a,
<a name="l00144"></a>00144                                 <a class="code" href="structfloat3.html">float3</a>* mfp_table_b,
<a name="l00145"></a>00145                                 <span class="keyword">struct</span> <a class="code" href="structrayleigh__struct.html" title="Structure storing the data of the Rayleigh interaction sampling model (equivalent to PENELOPE&#39;s commo...">rayleigh_struct</a>* rayleigh_table,
<a name="l00146"></a>00146                                 <span class="keyword">struct</span> <a class="code" href="structcompton__struct.html" title="Structure storing the data of the Compton interaction sampling model (equivalent to PENELOPE&#39;s common...">compton_struct</a>* compton_table,
<a name="l00147"></a>00147                                 <span class="keyword">struct</span> <a class="code" href="structdetector__struct.html" title="Structure storing the data defining the x-ray detector.">detector_struct</a>* detector_data_array,
<a name="l00148"></a>00148                                 <span class="keyword">struct</span> <a class="code" href="structsource__struct.html" title="Structure storing the data defining the source model (except for the energy spectrum).">source_struct</a>* source_data_array, 
<a name="l00149"></a>00149                                 <a class="code" href="structulonglong2.html">ulonglong2</a>* materials_dose)
<a name="l00150"></a>00150 <span class="preprocessor">#endif</span>
<a name="l00151"></a>00151 <span class="preprocessor"></span>{
<a name="l00152"></a>00152   <span class="comment">// -- Declare the track state variables:</span>
<a name="l00153"></a>00153   <a class="code" href="structfloat3.html">float3</a> position, direction;
<a name="l00154"></a>00154   <span class="keywordtype">float</span> energy, step, prob, randno, mfp_density, mfp_Woodcock;
<a name="l00155"></a>00155   <a class="code" href="structfloat3.html">float3</a> mfp_table_read_a, mfp_table_read_b;
<a name="l00156"></a>00156   <a class="code" href="structint2.html">int2</a> seed;
<a name="l00157"></a>00157   <span class="keywordtype">int</span> index;
<a name="l00158"></a>00158   <span class="keywordtype">int</span> material0,        <span class="comment">// Current material, starting at 0 for 1st material</span>
<a name="l00159"></a>00159       material_old;     <span class="comment">// Flag to mark a material or energy change</span>
<a name="l00160"></a>00160   <span class="keywordtype">signed</span> <span class="keywordtype">char</span> scatter_state;    <span class="comment">// Flag for scatter images: scatter_state=0 for non-scattered, =1 for Compton, =2 for Rayleigh, and =3 for multiple scatter.</span>
<a name="l00161"></a>00161 
<a name="l00162"></a>00162   <span class="comment">// -- Store the Compton table in shared memory from global memory:</span>
<a name="l00163"></a>00163   <span class="comment">//    For Compton and Rayleigh the access to memory is not coherent and the caching capability do not speeds up the accesses, they actually slows down the acces to other data.</span>
<a name="l00164"></a>00164 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00165"></a>00165 <span class="preprocessor"></span>  __shared__
<a name="l00166"></a>00166 <span class="preprocessor">#endif</span>
<a name="l00167"></a>00167 <span class="preprocessor"></span>  <span class="keyword">struct </span><a class="code" href="structcompton__struct.html" title="Structure storing the data of the Compton interaction sampling model (equivalent to PENELOPE&#39;s common...">compton_struct</a> cgco_SHARED;  
<a name="l00168"></a>00168 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00169"></a>00169 <span class="preprocessor"></span>  __shared__
<a name="l00170"></a>00170 <span class="preprocessor">#endif</span>
<a name="l00171"></a>00171 <span class="preprocessor"></span>  <span class="keyword">struct </span><a class="code" href="structdetector__struct.html" title="Structure storing the data defining the x-ray detector.">detector_struct</a> detector_data_SHARED;
<a name="l00172"></a>00172 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00173"></a>00173 <span class="preprocessor"></span>  __shared__
<a name="l00174"></a>00174 <span class="preprocessor">#endif </span>
<a name="l00175"></a>00175 <span class="preprocessor"></span>  <span class="keyword">struct </span><a class="code" href="structsource__struct.html" title="Structure storing the data defining the source model (except for the energy spectrum).">source_struct</a> source_data_SHARED;    
<a name="l00176"></a>00176 
<a name="l00177"></a>00177     
<a name="l00178"></a>00178 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00179"></a>00179 <span class="preprocessor"></span>  <span class="keywordflow">if</span> (0==threadIdx.x)  <span class="comment">// First GPU thread copies the variables to shared memory</span>
<a name="l00180"></a>00180   {
<a name="l00181"></a>00181 <span class="preprocessor">#endif</span>
<a name="l00182"></a>00182 <span class="preprocessor"></span>
<a name="l00183"></a>00183     <span class="comment">// -Copy the current source, detector data from global to shared memory for fast access:</span>
<a name="l00184"></a>00184     source_data_SHARED    = source_data_array[num_p];      
<a name="l00185"></a>00185     detector_data_SHARED  = detector_data_array[num_p];    <span class="comment">// Copy the long array to a single instance in shared memory for the current projection</span>
<a name="l00186"></a>00186         
<a name="l00187"></a>00187     <span class="comment">// -Copy the compton data to shared memory:</span>
<a name="l00188"></a>00188     cgco_SHARED = *compton_table;
<a name="l00189"></a>00189     
<a name="l00190"></a>00190 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00191"></a>00191 <span class="preprocessor"></span>  }
<a name="l00192"></a>00192   __syncthreads();     <span class="comment">// Make sure all threads will see the initialized shared variable  </span>
<a name="l00193"></a>00193 <span class="preprocessor">#endif</span>
<a name="l00194"></a>00194 <span class="preprocessor"></span>
<a name="l00195"></a>00195 
<a name="l00196"></a>00196   <span class="comment">// -- Initialize the RANECU generator in a position far away from the previous history:</span>
<a name="l00197"></a>00197 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00198"></a>00198 <span class="preprocessor"></span>  <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae7bd0b93b4d5eedb15822fa0ef622b48" title="Initialize the pseudo-random number generator (PRNG) RANECU to a position far away from the previous ...">init_PRNG</a>((threadIdx.x + blockIdx.x*blockDim.x), histories_per_thread, seed_input, &amp;seed);   <span class="comment">// Using a 1D block</span>
<a name="l00199"></a>00199 <span class="preprocessor">#else</span>
<a name="l00200"></a>00200 <span class="preprocessor"></span>  <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae7bd0b93b4d5eedb15822fa0ef622b48" title="Initialize the pseudo-random number generator (PRNG) RANECU to a position far away from the previous ...">init_PRNG</a>(history_batch, histories_per_thread, seed_input, &amp;seed);
<a name="l00201"></a>00201 <span class="preprocessor">#endif</span>
<a name="l00202"></a>00202 <span class="preprocessor"></span>
<a name="l00203"></a>00203   
<a name="l00204"></a>00204   <span class="comment">// -- Loop for the &quot;histories_per_thread&quot; particles in the current history_batch:</span>
<a name="l00205"></a>00205 
<a name="l00206"></a>00206   <span class="keywordflow">for</span>( ; histories_per_thread&gt;0; histories_per_thread--)
<a name="l00207"></a>00207   {
<a name="l00208"></a>00208         <span class="comment">//  printf(&quot;\n\n********* NEW HISTORY:  %d    [seeds: %d, %d]\n\n&quot;, histories_per_thread, seed.x, seed.y); //  fflush(stdout);  // !!Verbose!! calling printf from the GPU is possible but if multiple threads call it at the same time some output will be lost.</span>
<a name="l00209"></a>00209 
<a name="l00210"></a>00210     <span class="keywordtype">int</span> absvox = 1;
<a name="l00211"></a>00211     
<a name="l00212"></a>00212     <span class="comment">// -- Call the source function to get a primary x ray:</span>
<a name="l00213"></a>00213     <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a6807df36f6f2501a9f57ea4fce81b4f2" title="Source that creates primary x rays, according to the defined source model.">source</a>(&amp;position, &amp;direction, &amp;energy, &amp;seed, &amp;absvox, &amp;source_data_SHARED, &amp;detector_data_SHARED);
<a name="l00214"></a>00214 
<a name="l00215"></a>00215     scatter_state = (<span class="keywordtype">signed</span> char)0;     <span class="comment">// Reset previous scatter state: new non-scattered particle loaded</span>
<a name="l00216"></a>00216 
<a name="l00217"></a>00217     <span class="comment">// -- Find the current energy bin by truncation (this could be pre-calculated for a monoenergetic beam):    </span>
<a name="l00218"></a>00218     <span class="comment">//    The initialization host code made sure that the sampled energy will always be within the tabulated energies (index never negative or too large).</span>
<a name="l00219"></a>00219 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00220"></a>00220 <span class="preprocessor"></span>    index = __float2int_rd((energy-<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#ac07eade974afb1b75d16516d1eb55a33">e0</a>)*<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#a4b743ecf01ca7d075cc6b36b9ce2f0a1">ide</a>);  <span class="comment">// Using CUDA function to convert float to integer rounding down (towards minus infinite)</span>
<a name="l00221"></a>00221 <span class="preprocessor">#else</span>
<a name="l00222"></a>00222 <span class="preprocessor"></span>    index = (int)((energy-<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#ac07eade974afb1b75d16516d1eb55a33">e0</a>)*<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#a4b743ecf01ca7d075cc6b36b9ce2f0a1">ide</a> + 0.00001f);    <span class="comment">// Adding EPSILON to truncate to INT towards minus infinite. There may be a small error for energy&lt;=mfp_table_data_CONST.e0 but this case is irrelevant (particles will always have more energy than e0).</span>
<a name="l00223"></a>00223 <span class="preprocessor">#endif          </span>
<a name="l00224"></a>00224 <span class="preprocessor"></span>
<a name="l00225"></a>00225   
<a name="l00226"></a>00226     <span class="comment">// -- Get the minimum mfp at the current energy using linear interpolation (Woodcock tracking):      </span>
<a name="l00227"></a>00227     {
<a name="l00228"></a>00228       <a class="code" href="structfloat2.html">float2</a> mfp_Woodcock_read = mfp_Woodcock_table[index];   <span class="comment">// Read the 2 parameters for the linear interpolation in a single read from global memory</span>
<a name="l00229"></a>00229       mfp_Woodcock = mfp_Woodcock_read.<a class="code" href="structfloat2.html#a0bb68b4a8de04ffea5f7aae53c48a613">x</a> + energy * mfp_Woodcock_read.<a class="code" href="structfloat2.html#a0f1298f22c4ee20a369cd3b9c25b2cc6">y</a>;   <span class="comment">// Interpolated minimum MFP          </span>
<a name="l00230"></a>00230     }
<a name="l00231"></a>00231 
<a name="l00232"></a>00232 
<a name="l00233"></a>00233     <span class="comment">// -- Reset previous material to force a recalculation of the MFPs (negative materials are not allowed in the voxels):</span>
<a name="l00234"></a>00234     material_old  = -1;
<a name="l00235"></a>00235 
<a name="l00236"></a>00236     <span class="comment">// *** X-ray interaction loop:</span>
<a name="l00237"></a>00237     <span class="keywordflow">for</span>(;;)
<a name="l00238"></a>00238     {
<a name="l00239"></a>00239       
<a name="l00240"></a>00240       <span class="keywordflow">if</span> (absvox&lt;0)   <span class="comment">// !!DeBuG!!  MC-GPU_v1.3 ==&gt; if I move this &quot;if&quot; above the code runs much slower!? Why???</span>
<a name="l00241"></a>00241           <span class="keywordflow">break</span>;    <span class="comment">// -- Primary particle was not pointing to the voxel region! (but may still be detected after moving in vacuum in a straight line).      </span>
<a name="l00242"></a>00242 
<a name="l00243"></a>00243 
<a name="l00244"></a>00244       <span class="comment">// *** Virtual interaction loop:  // New loop structure in MC-GPU_v1.3: simulate all virtual events before sampling Compton &amp; Rayleigh:  // !!DeBuG!!</span>
<a name="l00245"></a>00245       
<a name="l00246"></a>00246       <a class="code" href="structfloat2.html">float2</a> matdens;
<a name="l00247"></a>00247       <a class="code" href="structshort3.html">short3</a> voxel_coord;    <span class="comment">// Variable used only by DOSE TALLY</span>
<a name="l00248"></a>00248 
<a name="l00249"></a>00249       <span class="keywordflow">do</span>
<a name="l00250"></a>00250       {     
<a name="l00251"></a>00251         step = -(mfp_Woodcock)*logf(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(&amp;seed));   <span class="comment">// Using the minimum MFP in the geometry for the input energy (Woodcock tracking)</span>
<a name="l00252"></a>00252           
<a name="l00253"></a>00253         position.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> += step*direction.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l00254"></a>00254         position.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> += step*direction.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l00255"></a>00255         position.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> += step*direction.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00256"></a>00256 
<a name="l00257"></a>00257         <span class="comment">// -- Locate the new particle in the voxel geometry:      </span>
<a name="l00258"></a>00258         absvox = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a562dba755f1abdb20b1088d3121f4326" title="Find the voxel that contains the current position.">locate_voxel</a>(&amp;position, &amp;voxel_coord);   <span class="comment">// Get the voxel number at the current position and the voxel coordinates (used to check if inside the dose ROI in DOSE TALLY).</span>
<a name="l00259"></a>00259         <span class="keywordflow">if</span> (absvox&lt;0)
<a name="l00260"></a>00260           <span class="keywordflow">break</span>;    <span class="comment">// -- Particle escaped the voxel region! (&quot;index&quot; is still &gt;0 at this moment)</span>
<a name="l00261"></a>00261           
<a name="l00262"></a>00262         matdens = voxel_mat_dens[absvox];     <span class="comment">// Get the voxel material and density in a single read from global memory</span>
<a name="l00263"></a>00263         material0 = (int)(matdens.<a class="code" href="structfloat2.html#a0bb68b4a8de04ffea5f7aae53c48a613">x</a> - 1);   <span class="comment">// Set the current material by truncation, and set 1st material to value &#39;0&#39;.</span>
<a name="l00264"></a>00264 
<a name="l00265"></a>00265         <span class="comment">// -- Get the data for the linear interpolation of the interaction MFPs, in case the energy or material have changed:</span>
<a name="l00266"></a>00266         <span class="keywordflow">if</span> (material0 != material_old)
<a name="l00267"></a>00267         {
<a name="l00268"></a>00268           mfp_table_read_a = mfp_table_a[index*(<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>)+material0];
<a name="l00269"></a>00269           mfp_table_read_b = mfp_table_b[index*(<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>)+material0];
<a name="l00270"></a>00270           material_old = material0;                                              <span class="comment">// Store the new material</span>
<a name="l00271"></a>00271         }
<a name="l00272"></a>00272         
<a name="l00273"></a>00273         <span class="comment">// *** Apply Woodcock tracking:</span>
<a name="l00274"></a>00274         mfp_density = mfp_Woodcock * matdens.<a class="code" href="structfloat2.html#a0f1298f22c4ee20a369cd3b9c25b2cc6">y</a>;
<a name="l00275"></a>00275         <span class="comment">// -- Calculate probability of delta scattering, using the total mean free path for the current material and energy (linear interpolation):</span>
<a name="l00276"></a>00276         prob = 1.0f - mfp_density * (mfp_table_read_a.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + energy * mfp_table_read_b.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>);
<a name="l00277"></a>00277         randno = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(&amp;seed);    <span class="comment">// Sample uniform PRN</span>
<a name="l00278"></a>00278       }
<a name="l00279"></a>00279       <span class="keywordflow">while</span> (randno&lt;prob);   <span class="comment">// [Iterate if there is a delta scattering event]</span>
<a name="l00280"></a>00280 
<a name="l00281"></a>00281       <span class="keywordflow">if</span> (absvox&lt;0)
<a name="l00282"></a>00282         <span class="keywordflow">break</span>;    <span class="comment">// -- Particle escaped the voxel region! Break the interaction loop to call tally image.</span>
<a name="l00283"></a>00283 
<a name="l00284"></a>00284         
<a name="l00285"></a>00285       <span class="comment">// The GPU threads will be stopped and waiting here until ALL threads have a REAL event: </span>
<a name="l00286"></a>00286 
<a name="l00287"></a>00287       <span class="comment">// -- Real event takes place! Check the kind of event and sample the effects of the interaction:</span>
<a name="l00288"></a>00288       
<a name="l00289"></a>00289       prob += mfp_density * (mfp_table_read_a.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> + energy * mfp_table_read_b.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>);    <span class="comment">// Interpolate total Compton MFP (&#39;y&#39; component)</span>
<a name="l00290"></a>00290       <span class="keywordflow">if</span> (randno&lt;prob)   <span class="comment">// [Checking Compton scattering]</span>
<a name="l00291"></a>00291       {
<a name="l00292"></a>00292         <span class="comment">// *** Compton interaction:</span>
<a name="l00293"></a>00293 
<a name="l00294"></a>00294         <span class="comment">//  -- Sample new direction and energy:</span>
<a name="l00295"></a>00295         <span class="keywordtype">double</span> costh_Compton;
<a name="l00296"></a>00296         randno = energy;     <span class="comment">// Save temporal copy of the particle energy (variable randno not necessary until next sampling). DOSE TALLY</span>
<a name="l00297"></a>00297         
<a name="l00298"></a>00298         <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a6f212862e7fba1553ec9495f05443b2c" title="Random sampling of incoherent (Compton) scattering of photons, using the sampling algorithm from PENE...">GCOa</a>(&amp;energy, &amp;costh_Compton, &amp;material0, &amp;seed, &amp;cgco_SHARED);
<a name="l00299"></a>00299         <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a94ad9ed794eb8fd696e7e357ae391765" title="Rotates a vector; the rotation is specified by giving the polar and azimuthal angles in the &quot;self-fra...">rotate_double</a>(&amp;direction, costh_Compton, <span class="comment">/*phi=2*pi*PRN=*/</span> 6.28318530717958647693*<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(&amp;seed));
<a name="l00300"></a>00300 
<a name="l00301"></a>00301         randno = energy - randno;   <span class="comment">// Save temporal copy of the negative of the energy lost in the interaction.  DOSE TALLY</span>
<a name="l00302"></a>00302 
<a name="l00303"></a>00303         <span class="comment">// -- Find the new energy interval:</span>
<a name="l00304"></a>00304 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00305"></a>00305 <span class="preprocessor"></span>        index = __float2int_rd((energy-<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#ac07eade974afb1b75d16516d1eb55a33">e0</a>)*<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#a4b743ecf01ca7d075cc6b36b9ce2f0a1">ide</a>);  <span class="comment">// Using CUDA function to convert float to integer rounding down (towards minus infinite)</span>
<a name="l00306"></a>00306 <span class="preprocessor">#else</span>
<a name="l00307"></a>00307 <span class="preprocessor"></span>        index = (int)((energy-<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#ac07eade974afb1b75d16516d1eb55a33">e0</a>)*<a class="code" href="MC-GPU__v1_83_8h.html#a8a3e410fb23d04f8a27d3d71c2381d24" title="Global variable to be stored in the GPU constant memory defining the linear interpolation data...">mfp_table_data_CONST</a>.<a class="code" href="structlinear__interp.html#a4b743ecf01ca7d075cc6b36b9ce2f0a1">ide</a> + 0.00001f);    <span class="comment">// Adding EPSILON to truncate to INT</span>
<a name="l00308"></a>00308 <span class="preprocessor">#endif          </span>
<a name="l00309"></a>00309 <span class="preprocessor"></span>
<a name="l00310"></a>00310         
<a name="l00311"></a>00311         <span class="keywordflow">if</span> (index&gt;-1)  <span class="comment">// &#39;index&#39; will be negative only when the energy is below the tabulated minimum energy: particle will be then absorbed (rejected) after tallying the dose.</span>
<a name="l00312"></a>00312         {          
<a name="l00313"></a>00313           <span class="comment">// -- Get the Woodcock MFP for the new energy (energy above minimum cutoff):</span>
<a name="l00314"></a>00314           <a class="code" href="structfloat2.html">float2</a> mfp_Woodcock_read = mfp_Woodcock_table[index];   <span class="comment">// Read the 2 parameters for the linear interpolation in a single read from global memory</span>
<a name="l00315"></a>00315           mfp_Woodcock = mfp_Woodcock_read.<a class="code" href="structfloat2.html#a0bb68b4a8de04ffea5f7aae53c48a613">x</a> + energy * mfp_Woodcock_read.<a class="code" href="structfloat2.html#a0f1298f22c4ee20a369cd3b9c25b2cc6">y</a>;   <span class="comment">// Interpolated minimum MFP</span>
<a name="l00316"></a>00316 
<a name="l00317"></a>00317           material_old = -2;    <span class="comment">// Set an impossible material to force an update of the MFPs data for the nex energy interval</span>
<a name="l00318"></a>00318 
<a name="l00319"></a>00319           <span class="comment">// -- Update scatter state:</span>
<a name="l00320"></a>00320           <span class="keywordflow">if</span> (scatter_state==(<span class="keywordtype">signed</span> <span class="keywordtype">char</span>)0)
<a name="l00321"></a>00321             scatter_state = (<span class="keywordtype">signed</span> <span class="keywordtype">char</span>)1;   <span class="comment">// Set scatter_state == 1: Compton scattered particle</span>
<a name="l00322"></a>00322           <span class="keywordflow">else</span>
<a name="l00323"></a>00323             scatter_state = (<span class="keywordtype">signed</span> char)3;   <span class="comment">// Set scatter_state == 3: Multi-scattered particle</span>
<a name="l00324"></a>00324         }
<a name="l00325"></a>00325 
<a name="l00326"></a>00326       }
<a name="l00327"></a>00327       <span class="keywordflow">else</span>
<a name="l00328"></a>00328       {
<a name="l00329"></a>00329         prob += mfp_density * (mfp_table_read_a.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> + energy * mfp_table_read_b.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>);    <span class="comment">// Interpolate total Rayleigh MFP (&#39;z&#39; component)</span>
<a name="l00330"></a>00330         <span class="keywordflow">if</span> (randno&lt;prob)   <span class="comment">// [Checking Rayleigh scattering]</span>
<a name="l00331"></a>00331         {
<a name="l00332"></a>00332           <span class="comment">// *** Rayleigh interaction:</span>
<a name="l00333"></a>00333 
<a name="l00334"></a>00334           <span class="comment">//  -- Sample angular deflection:</span>
<a name="l00335"></a>00335           <span class="keywordtype">double</span> costh_Rayleigh;
<a name="l00336"></a>00336           <span class="keywordtype">float</span> pmax_current = rayleigh_table-&gt;<a class="code" href="structrayleigh__struct.html#a3030ed3f1f07daf452bf53fe6aa8f7a1">pmax</a>[(index+1)*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>+material0];   <span class="comment">// Get max (ie, value for next bin?) cumul prob square form factor for Rayleigh sampling</span>
<a name="l00337"></a>00337 
<a name="l00338"></a>00338           <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a5cafd4ebb3d155746b16ca7d231ecad1" title="Sample a Rayleigh interaction using the sampling algorithm used in PENELOPE 2006.">GRAa</a>(&amp;energy, &amp;costh_Rayleigh, &amp;material0, &amp;pmax_current, &amp;seed, rayleigh_table);
<a name="l00339"></a>00339           <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a94ad9ed794eb8fd696e7e357ae391765" title="Rotates a vector; the rotation is specified by giving the polar and azimuthal angles in the &quot;self-fra...">rotate_double</a>(&amp;direction, costh_Rayleigh, <span class="comment">/*phi=2*pi*PRN=*/</span> 6.28318530717958647693*<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(&amp;seed));
<a name="l00340"></a>00340 
<a name="l00341"></a>00341           <span class="comment">// -- Update scatter state:</span>
<a name="l00342"></a>00342           <span class="keywordflow">if</span> (scatter_state==(<span class="keywordtype">signed</span> <span class="keywordtype">char</span>)0)
<a name="l00343"></a>00343             scatter_state = (<span class="keywordtype">signed</span> char)2;   <span class="comment">// Set scatter_state == 1: Rayleigh scattered particle</span>
<a name="l00344"></a>00344           <span class="keywordflow">else</span>
<a name="l00345"></a>00345             scatter_state = (<span class="keywordtype">signed</span> char)3;   <span class="comment">// Set scatter_state == 3: Multi-scattered particle</span>
<a name="l00346"></a>00346 
<a name="l00347"></a>00347         }
<a name="l00348"></a>00348         <span class="keywordflow">else</span>
<a name="l00349"></a>00349         {
<a name="l00350"></a>00350           <span class="comment">// *** Photoelectric interaction (or pair production): mark particle for absorption after dose tally (ie, index&lt;0)!</span>
<a name="l00351"></a>00351           randno = -energy;   <span class="comment">// Save temporal copy of the (negative) energy deposited in the interaction (variable randno not necessary anymore).</span>
<a name="l00352"></a>00352           index = -11;       <span class="comment">// A negative &quot;index&quot; marks that the particle was absorved and that it will never arrive at the detector.</span>
<a name="l00353"></a>00353         }
<a name="l00354"></a>00354       }
<a name="l00355"></a>00355     
<a name="l00356"></a>00356       <span class="comment">//  -- Tally the dose deposited in Compton and photoelectric interactions:</span>
<a name="l00357"></a>00357       <span class="keywordflow">if</span> (randno&lt;-0.001f)
<a name="l00358"></a>00358       {
<a name="l00359"></a>00359         <span class="keywordtype">float</span> Edep = -1.0f*randno;   <span class="comment">// If any energy was deposited, this variable will temporarily store the negative value of Edep.</span>
<a name="l00360"></a>00360         
<a name="l00361"></a>00361         <span class="comment">//  -- Tally the dose deposited in the current material, if enabled (ie, array allocated and not null):</span>
<a name="l00362"></a>00362         <span class="keywordflow">if</span> (materials_dose!=NULL)
<a name="l00363"></a>00363           <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a45a0719009ba55cca7af72f08f1090fb" title="Tally the depose deposited inside each material.">tally_materials_dose</a>(&amp;Edep, &amp;material0, materials_dose);    <span class="comment">// !!tally_materials_dose!!</span>
<a name="l00364"></a>00364 
<a name="l00365"></a>00365         <span class="comment">//  -- Tally the energy deposited in the current voxel, if enabled (tally disabled when dose_ROI_x_max_CONST is negative). DOSE TALLY</span>
<a name="l00366"></a>00366         <span class="keywordflow">if</span> (<a class="code" href="MC-GPU__v1_83_8h.html#a3c22177ce0c73c17ca8716d997678b61">dose_ROI_x_max_CONST</a> &gt; -1)
<a name="l00367"></a>00367           <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aac6cb77943874a858ac5a79cbd58f1e7" title="Tally the dose deposited in the voxels.">tally_voxel_energy_deposition</a>(&amp;Edep, &amp;voxel_coord, voxels_Edep);
<a name="l00368"></a>00368 
<a name="l00369"></a>00369       }    
<a name="l00370"></a>00370 
<a name="l00371"></a>00371       <span class="comment">// -- Break interaction loop for particles that have been absorved or with energy below the tabulated cutoff: particle is &quot;absorbed&quot; (ie, track discontinued).</span>
<a name="l00372"></a>00372       <span class="keywordflow">if</span> (index&lt;0)
<a name="l00373"></a>00373         <span class="keywordflow">break</span>;  
<a name="l00374"></a>00374       
<a name="l00375"></a>00375     }   <span class="comment">// [Cycle the X-ray interaction loop]</span>
<a name="l00376"></a>00376 
<a name="l00377"></a>00377     <span class="keywordflow">if</span> (index&gt;-1)
<a name="l00378"></a>00378     {
<a name="l00379"></a>00379       <span class="comment">// -- Particle escaped the voxels but was not absorbed, check if it will arrive at the detector and tally its energy:</span>
<a name="l00380"></a>00380       <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aac14e3b20e073c6b169f2bd5c23da281" title="Tally a radiographic projection image.">tally_image</a>(&amp;energy, &amp;position, &amp;direction, &amp;scatter_state, image, &amp;source_data_SHARED, &amp;detector_data_SHARED);
<a name="l00381"></a>00381     }
<a name="l00382"></a>00382   }   <span class="comment">// [Continue with a new history]</span>
<a name="l00383"></a>00383 
<a name="l00384"></a>00384 }   <span class="comment">// [All tracks simulated for this kernel call: return to CPU]</span>
<a name="l00385"></a>00385 
<a name="l00386"></a>00386 
<a name="l00387"></a>00387 
<a name="l00388"></a>00388 
<a name="l00389"></a>00389 
<a name="l00390"></a>00390 <span class="comment"></span>
<a name="l00391"></a>00391 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00392"></a>00392 <span class="comment">//!  Tally the dose deposited in the voxels.</span>
<a name="l00393"></a>00393 <span class="comment">//!  This function is called whenever a particle suffers a Compton or photoelectric</span>
<a name="l00394"></a>00394 <span class="comment">//!  interaction. It is not necessary to call this function if the dose tally</span>
<a name="l00395"></a>00395 <span class="comment">//!  was disabled in the input file (ie, dose_ROI_x_max_CONST &lt; 0).</span>
<a name="l00396"></a>00396 <span class="comment">//!  Electrons are not transported in MC-GPU and therefore we are approximating</span>
<a name="l00397"></a>00397 <span class="comment">//!  that the dose is equal to the KERMA (energy released by the photons alone).</span>
<a name="l00398"></a>00398 <span class="comment">//!  This approximation is acceptable when there is electronic equilibrium and when</span>
<a name="l00399"></a>00399 <span class="comment">//!  the range of the secondary electrons is shorter than the voxel size. Usually the</span>
<a name="l00400"></a>00400 <span class="comment">//!  doses will be acceptable for photon energies below 1 MeV. The dose estimates may</span>
<a name="l00401"></a>00401 <span class="comment">//!  not be accurate at the interface of low density volumes.</span>
<a name="l00402"></a>00402 <span class="comment">//!</span>
<a name="l00403"></a>00403 <span class="comment">//!  We need to use atomicAdd() in the GPU to prevent that multiple threads update the </span>
<a name="l00404"></a>00404 <span class="comment">//!  same voxel at the same time, which would result in a lose of information.</span>
<a name="l00405"></a>00405 <span class="comment">//!  This is very improbable when using a large number of voxels but gives troubles </span>
<a name="l00406"></a>00406 <span class="comment">//!  with a simple geometries with few voxels (in this case the atomicAdd will slow </span>
<a name="l00407"></a>00407 <span class="comment">//!  down the code because threads will update the voxel dose secuentially).</span>
<a name="l00408"></a>00408 <span class="comment">//!</span>
<a name="l00409"></a>00409 <span class="comment">//!</span>
<a name="l00410"></a>00410 <span class="comment">//!       @param[in] Edep   Energy deposited in the interaction</span>
<a name="l00411"></a>00411 <span class="comment">//!       @param[in] voxel_coord   Voxel coordinates, needed to check if particle located inside the input region of interest (ROI)</span>
<a name="l00412"></a>00412 <span class="comment">//!       @param[out] voxels_Edep   ulonglong2 array containing the 3D voxel dose and dose^2 (ie, uncertainty) as unsigned integers scaled by SCALE_eV.</span>
<a name="l00413"></a>00413 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00414"></a>00414 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00415"></a>00415 <span class="preprocessor"></span>__device__
<a name="l00416"></a>00416 <span class="preprocessor">#endif</span>
<a name="l00417"></a>00417 <span class="preprocessor"></span><span class="keyword">inline</span> 
<a name="l00418"></a><a class="code" href="MC-GPU__v1_83_8h.html#abf1491ec13d6068d2753516a3721cfe2">00418</a> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aac6cb77943874a858ac5a79cbd58f1e7" title="Tally the dose deposited in the voxels.">tally_voxel_energy_deposition</a>(<span class="keywordtype">float</span>* Edep, <a class="code" href="structshort3.html">short3</a>* voxel_coord, <a class="code" href="structulonglong2.html">ulonglong2</a>* voxels_Edep)
<a name="l00419"></a>00419 {
<a name="l00420"></a>00420 
<a name="l00421"></a>00421     <span class="comment">// !!DeBuG!! Maybe it would be faster to store a 6 element struct and save temp copy?? struct_short_int_x6_align16  dose_ROI_size = dose_ROI_size_CONST;   // Get ROI coordinates from GPU constant memory and store temporal copy</span>
<a name="l00422"></a>00422   
<a name="l00423"></a>00423   <span class="keywordflow">if</span>((voxel_coord-&gt;<a class="code" href="structshort3.html#a5d5d88cacca4de50a2334be0033a126a">x</a> &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a5d34c13b457e76dec83755628ead56fc" title="Global variable to be stored in the GPU constant memory defining the coordinates of the dose depositi...">dose_ROI_x_min_CONST</a>) || (voxel_coord-&gt;<a class="code" href="structshort3.html#a5d5d88cacca4de50a2334be0033a126a">x</a> &gt; <a class="code" href="MC-GPU__v1_83_8h.html#a3c22177ce0c73c17ca8716d997678b61">dose_ROI_x_max_CONST</a>) ||
<a name="l00424"></a>00424      (voxel_coord-&gt;<a class="code" href="structshort3.html#a6056b3f7d064291aaf9121575ec83866">y</a> &lt; <a class="code" href="MC-GPU__v1_83_8h.html#ac161c8f74fd460e82e73917b3afff5fc">dose_ROI_y_min_CONST</a>) || (voxel_coord-&gt;<a class="code" href="structshort3.html#a6056b3f7d064291aaf9121575ec83866">y</a> &gt; <a class="code" href="MC-GPU__v1_83_8h.html#ab591ece7d8e024dcc8de71f59bc1cea5">dose_ROI_y_max_CONST</a>) ||
<a name="l00425"></a>00425      (voxel_coord-&gt;<a class="code" href="structshort3.html#a803658b690ee5726c60a2b99b99b9143">z</a> &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a81db18966d58018976252f898cd905d9">dose_ROI_z_min_CONST</a>) || (voxel_coord-&gt;<a class="code" href="structshort3.html#a803658b690ee5726c60a2b99b99b9143">z</a> &gt; <a class="code" href="MC-GPU__v1_83_8h.html#aea8fcb598ccc9c18ae601f9c7d936b9d">dose_ROI_z_max_CONST</a>))
<a name="l00426"></a>00426     {
<a name="l00427"></a>00427       <span class="keywordflow">return</span>;   <span class="comment">// -- Particle outside the ROI: return without tallying anything.</span>
<a name="l00428"></a>00428     }
<a name="l00429"></a>00429 
<a name="l00430"></a>00430   <span class="comment">// -- Particle inside the ROI: tally Edep.</span>
<a name="l00431"></a>00431   <span class="keyword">register</span> <span class="keywordtype">int</span> DX = 1 + (int)(<a class="code" href="MC-GPU__v1_83_8h.html#a3c22177ce0c73c17ca8716d997678b61">dose_ROI_x_max_CONST</a> - <a class="code" href="MC-GPU__v1_83_8h.html#a5d34c13b457e76dec83755628ead56fc" title="Global variable to be stored in the GPU constant memory defining the coordinates of the dose depositi...">dose_ROI_x_min_CONST</a>);
<a name="l00432"></a>00432   <span class="keyword">register</span> <span class="keywordtype">int</span> num_voxel = (int)(voxel_coord-&gt;<a class="code" href="structshort3.html#a5d5d88cacca4de50a2334be0033a126a">x</a>-<a class="code" href="MC-GPU__v1_83_8h.html#a5d34c13b457e76dec83755628ead56fc" title="Global variable to be stored in the GPU constant memory defining the coordinates of the dose depositi...">dose_ROI_x_min_CONST</a>) + ((int)(voxel_coord-&gt;<a class="code" href="structshort3.html#a6056b3f7d064291aaf9121575ec83866">y</a>-<a class="code" href="MC-GPU__v1_83_8h.html#ac161c8f74fd460e82e73917b3afff5fc">dose_ROI_y_min_CONST</a>))*DX + ((<span class="keywordtype">int</span>)(voxel_coord-&gt;<a class="code" href="structshort3.html#a803658b690ee5726c60a2b99b99b9143">z</a>-<a class="code" href="MC-GPU__v1_83_8h.html#a81db18966d58018976252f898cd905d9">dose_ROI_z_min_CONST</a>))*DX*(1 + (int)(<a class="code" href="MC-GPU__v1_83_8h.html#ab591ece7d8e024dcc8de71f59bc1cea5">dose_ROI_y_max_CONST</a>-<a class="code" href="MC-GPU__v1_83_8h.html#ac161c8f74fd460e82e73917b3afff5fc">dose_ROI_y_min_CONST</a>));
<a name="l00433"></a>00433   
<a name="l00434"></a>00434 <span class="preprocessor">   #ifdef USING_CUDA</span>
<a name="l00435"></a>00435 <span class="preprocessor"></span>     atomicAdd(&amp;voxels_Edep[num_voxel].x, __float2ull_rn((*Edep)*<a class="code" href="MC-GPU__v1_83_8h.html#ac0c03cf5e803404f5dc2a42467395a3a" title="Value to scale the deposited energy in the pixels so that it can be stored as a long long integer ins...">SCALE_eV</a>) );    <span class="comment">// Energy deposited at the voxel, scaled by the factor SCALE_eV and rounded.</span>
<a name="l00436"></a>00436      atomicAdd(&amp;voxels_Edep[num_voxel].y, __float2ull_rn((*Edep)*(*Edep)) );     <span class="comment">// (not using SCALE_eV for std_dev to prevent overflow)           </span>
<a name="l00437"></a>00437 <span class="preprocessor">   #else</span>
<a name="l00438"></a>00438 <span class="preprocessor"></span>     voxels_Edep[num_voxel].<a class="code" href="structulonglong2.html#a0d94e96b26765bf28a1a4b76391b88bb">x</a> += (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)((*Edep)*SCALE_eV + 0.5f);
<a name="l00439"></a>00439      voxels_Edep[num_voxel].<a class="code" href="structulonglong2.html#a8453eaa1efda22c449296aa2ebcf0c7f">y</a> += (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)((*Edep)*(*Edep) + 0.5f);
<a name="l00440"></a>00440 <span class="preprocessor">   #endif</span>
<a name="l00441"></a>00441 <span class="preprocessor"></span>          
<a name="l00442"></a>00442   <span class="keywordflow">return</span>;
<a name="l00443"></a>00443 }
<a name="l00444"></a>00444 
<a name="l00445"></a>00445 <span class="comment"></span>
<a name="l00446"></a>00446 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00447"></a>00447 <span class="comment">//!  Tally a radiographic projection image.</span>
<a name="l00448"></a>00448 <span class="comment">//!  This function is called whenever a particle escapes the voxelized volume.</span>
<a name="l00449"></a>00449 <span class="comment">//!  The code checks if the particle would arrive at the detector if it kept</span>
<a name="l00450"></a>00450 <span class="comment">//!  moving in a straight line after exiting the voxels (assuming vacuum enclosure).</span>
<a name="l00451"></a>00451 <span class="comment">//!  An ideal image formation model is implemented: each pixel counts the total energy</span>
<a name="l00452"></a>00452 <span class="comment">//!  of the x rays that enter the pixel (100% detection efficiency for any energy).</span>
<a name="l00453"></a>00453 <span class="comment">//!  The image due to primaries and different kinds of scatter is tallied separately.</span>
<a name="l00454"></a>00454 <span class="comment">//!</span>
<a name="l00455"></a>00455 <span class="comment">//!  In the GPU, and atomicAdd() function is used to make sure that multiple threads do</span>
<a name="l00456"></a>00456 <span class="comment">//!  not update the same pixel at the same time, which would result in a lose of information.</span>
<a name="l00457"></a>00457 <span class="comment">//!  Since the atomicAdd function is only available for &#39;unsigned long long int&#39; data,</span>
<a name="l00458"></a>00458 <span class="comment">//!  the float pixel values are scaled by a factor &quot;SCALE_eV&quot; defined in the header file</span>
<a name="l00459"></a>00459 <span class="comment">//!  (eg, #define SCALE_eV 10000.0f) and stored as unsigned long long integers in main</span>
<a name="l00460"></a>00460 <span class="comment">//!  memory.</span>
<a name="l00461"></a>00461 <span class="comment">//!</span>
<a name="l00462"></a>00462 <span class="comment">//!  WARNING! If the total tallied signal (for all particles) is larger than &quot;1.8e19/SCALE_eV&quot;,</span>
<a name="l00463"></a>00463 <span class="comment">//!    there will be a bit overflow and the value will be reset to 0 giving bogus results.</span>
<a name="l00464"></a>00464 <span class="comment">//!</span>
<a name="l00465"></a>00465 <span class="comment">//!  WARNING! The detector plane should be located outside the voxels bounding box. However, since</span>
<a name="l00466"></a>00466 <span class="comment">//!    the particles are moved outside the bbox in the last step, they could cross the detector </span>
<a name="l00467"></a>00467 <span class="comment">//!    plane anyway. If the particles are less than 2.0 cm behind the detector, they are moved </span>
<a name="l00468"></a>00468 <span class="comment">//!    back and detected. Therefore the detector can be a few cm inside the bbox and still work.</span>
<a name="l00469"></a>00469 <span class="comment">//!    If the Woodcock mean free path is larger than the distance from the bbox to the detector, </span>
<a name="l00470"></a>00470 <span class="comment">//!    we may lose some particles behind the detector!</span>
<a name="l00471"></a>00471 <span class="comment">//!</span>
<a name="l00472"></a>00472 <span class="comment">//!</span>
<a name="l00473"></a>00473 <span class="comment">//!       @param[in] energy   X-ray energy</span>
<a name="l00474"></a>00474 <span class="comment">//!       @param[in] position   Particle position</span>
<a name="l00475"></a>00475 <span class="comment">//!       @param[in] direction   Particle direction (cosine vectors)</span>
<a name="l00476"></a>00476 <span class="comment">//!       @param[in] scatter_state  Flag marking primaries, single Compton, single Rayleigh or multiple scattered radiation</span>
<a name="l00477"></a>00477 <span class="comment">//!       @param[out] image   Integer array containing the image, ie, the pixel values (in tenths of meV)</span>
<a name="l00478"></a>00478 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00479"></a>00479 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00480"></a>00480 <span class="preprocessor"></span>__device__
<a name="l00481"></a>00481 <span class="preprocessor">#endif</span>
<a name="l00482"></a><a class="code" href="MC-GPU__v1_83_8h.html#aac14e3b20e073c6b169f2bd5c23da281">00482</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aac14e3b20e073c6b169f2bd5c23da281" title="Tally a radiographic projection image.">tally_image</a>(<span class="keywordtype">float</span>* energy, <a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>, <a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>, <span class="keywordtype">signed</span> <span class="keywordtype">char</span>* scatter_state, <span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span>* image, <span class="keyword">struct</span> <a class="code" href="structsource__struct.html" title="Structure storing the data defining the source model (except for the energy spectrum).">source_struct</a>* source_data_SHARED, <span class="keyword">struct</span> <a class="code" href="structdetector__struct.html" title="Structure storing the data defining the x-ray detector.">detector_struct</a>* detector_data_SHARED)
<a name="l00483"></a>00483 {
<a name="l00484"></a>00484   <span class="keywordtype">float</span> dist_detector, rotated_position;
<a name="l00485"></a>00485 
<a name="l00486"></a>00486   <span class="keywordflow">if</span> (detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#adbe41a01886ae097dac0fa908551d035">rotation_flag</a> == 1)    <span class="comment">// --&gt;  Initial source direction is not (0,1,0): detector has to be rotated to +Y to find the pixel number</span>
<a name="l00487"></a>00487   {
<a name="l00488"></a>00488     
<a name="l00489"></a>00489     <span class="comment">// *** Skip particles not moving towards the detector. </span>
<a name="l00490"></a>00490     <span class="comment">//       (a) Skip particles that were deflected more than 90 deg from the original source direction (backscatter).</span>
<a name="l00491"></a>00491     <span class="comment">//       (b) Skip particles located more than 10 cm behind the detector </span>
<a name="l00492"></a>00492     <span class="comment">//       (c) Skip particles for which the direction to the detector is way bigger than SDD (likely to intersect the plane outside the pixel region).</span>
<a name="l00493"></a>00493                   <span class="comment">// !!DeBuG!! NOTE: This may give problems for big detectors very close to the source</span>
<a name="l00494"></a>00494                   
<a name="l00495"></a>00495     <span class="comment">//      !!DeBuG!! Particles located after the detector will be moved back to the surface of the detector, but 10 cm maximum!!</span>
<a name="l00496"></a>00496     <span class="comment">//                In this way the detector can intersect the voxels bbox or be located right on the surface of the bbox: the particles will be </span>
<a name="l00497"></a>00497     <span class="comment">//                transported across the detector and until a little after the end of the bbox in the last step, but then moved back.</span>
<a name="l00498"></a>00498     <span class="comment">//                This algorithm will give correct results ONLY when the detector intersects just slightly the air space around the phantom,</span>
<a name="l00499"></a>00499     <span class="comment">//                so that the interactions after the detector are not significant (this happens sometimes using oblique beams).</span>
<a name="l00500"></a>00500     <span class="comment">//                I could remove particles after the detector using &quot;if (dist_detector&lt;0.0f) return;&quot;.</span>
<a name="l00501"></a>00501 
<a name="l00502"></a>00502     <span class="comment">//  (a) Calculate the angle between the particle and the initial direction (dot product): reject particle if cos_angle &lt; cos(89)==0 (angle&gt;89deg):</span>
<a name="l00503"></a>00503     <span class="comment">//      [Extra parenthesis are coded to suggest to the compiler the use of intrinsic multiply-add operations].</span>
<a name="l00504"></a>00504 
<a name="l00505"></a>00505     <span class="keyword">register</span> <span class="keywordtype">float</span> cos_angle = direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> * source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> +
<a name="l00506"></a>00506                               (direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> * source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> +
<a name="l00507"></a>00507                               (direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> * source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>));    
<a name="l00508"></a>00508     <span class="keywordflow">if</span> (cos_angle &lt; 0.025f)
<a name="l00509"></a>00509       <span class="keywordflow">return</span>;  <span class="comment">// Reject particle: Angle larger than 89 deg --&gt; particle moving parallel to the detector or backwards towards the source!</span>
<a name="l00510"></a>00510 
<a name="l00511"></a>00511     <span class="comment">//   (b) Find the distance from the current particle location (likely just after the surface of the voxel bbox) to the intersection with the detector plane:</span>
<a name="l00512"></a>00512     dist_detector = ( source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> * (detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1f0ce65af62c45925dfc96de3b796474">center</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> - position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>) +
<a name="l00513"></a>00513                      (source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> * (detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1f0ce65af62c45925dfc96de3b796474">center</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> - position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>) +
<a name="l00514"></a>00514                      (source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> * (detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1f0ce65af62c45925dfc96de3b796474">center</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> - position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>))) ) / cos_angle;
<a name="l00515"></a>00515 
<a name="l00516"></a>00516                         
<a name="l00517"></a>00517                      
<a name="l00518"></a>00518 <span class="comment">// !!DeBuG!!  IF&#39;s below (used in v1.2) are not needed when checking the x ray angle:</span>
<a name="l00519"></a>00519 <span class="comment">//   if (dist_detector &lt; -10.0f)   // !!DeBuG!! Is 10 cm enough or too much? Should I use 0? or allow any distance?</span>
<a name="l00520"></a>00520 <span class="comment">//      return;  // !!DeBuG!! Reject particles located more than 10 cm behind the detector. 10 cm was selected arbitrarily. Woodcock MFP for x-rays in bone: MFP 200 keV photons in bone ==&gt; 4 cm.</span>
<a name="l00521"></a>00521 <span class="comment">//      </span>
<a name="l00522"></a>00522 <span class="comment">//    if (fabsf(dist_detector)&gt;(2.1f*detector_data_CONST.sdd))          </span>
<a name="l00523"></a>00523 <span class="comment">//      return;  // Reject particle: distance to the detector plane too large, the particle is likely to travel almost parallel to the detector and will not be detected.</span>
<a name="l00524"></a>00524 
<a name="l00525"></a>00525             
<a name="l00526"></a>00526     <span class="comment">// *** Translate the particle to the detector plane (we assume the detector is completely absorbent: 100% detection efficiency):</span>
<a name="l00527"></a>00527     position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + dist_detector * direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l00528"></a>00528     position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> + dist_detector * direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l00529"></a>00529     position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> + dist_detector * direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00530"></a>00530 
<a name="l00531"></a>00531     <span class="comment">// *** Rotate the particle position vector to the default reference system where the detector is perpendicular to the +Y axis, then find out if the particle is located inside a pixel:</span>
<a name="l00532"></a>00532 <span class="preprocessor">    #ifdef USING_CUDA</span>
<a name="l00533"></a>00533 <span class="preprocessor"></span>      rotated_position = detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[0]*position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[1]*position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[2]*position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;  <span class="comment">// X coordinate</span>
<a name="l00534"></a>00534       <span class="keywordtype">int</span> pixel_coord_x = __float2int_rd((rotated_position - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>) * detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a3c6f9e59e89db06211ccea63287e2826">inv_pixel_size_X</a>);    <span class="comment">// Using CUDA intrinsic function to convert float to integer rounding down (towards minus infinite)</span>
<a name="l00535"></a>00535       <span class="keywordflow">if</span> ((pixel_coord_x&gt;-1)&amp;&amp;(pixel_coord_x&lt;detector_data_SHARED-&gt;num_pixels.x))
<a name="l00536"></a>00536       {
<a name="l00537"></a>00537         rotated_position = detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[6]*position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[7]*position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[8]*position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;  <span class="comment">// Z coordinate</span>
<a name="l00538"></a>00538         <span class="keywordtype">int</span> pixel_coord_z = __float2int_rd((rotated_position - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>) * detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#aa91fc94f674ad7af8eca577f4d8cd22b">inv_pixel_size_Z</a>);
<a name="l00539"></a>00539         <span class="keywordflow">if</span> ((pixel_coord_z&gt;-1)&amp;&amp;(pixel_coord_z&lt;detector_data_SHARED-&gt;num_pixels.y))
<a name="l00540"></a>00540         {
<a name="l00541"></a>00541           <span class="comment">// -- Particle enters the detector! Tally the particle energy in the corresponding pixel (in tenths of meV):</span>
<a name="l00542"></a>00542           <span class="comment">//    Using a CUDA atomic function (not available for global floats yet) to read and increase the pixel value in a single instruction, blocking interferences from other threads.</span>
<a name="l00543"></a>00543           <span class="comment">//    The offset for the primaries or scatter images are calculated considering that:</span>
<a name="l00544"></a>00544           <span class="comment">//      scatter_state=0 for non-scattered, =1 for Compton, =2 for Rayleigh, and =3 for multiple scatter.</span>
<a name="l00545"></a>00545           atomicAdd(( image +                                                               <span class="comment">// Pointer to beginning of image array</span>
<a name="l00546"></a>00546                     (<span class="keywordtype">int</span>)(*scatter_state) * detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#ac26d0b66ebf521346d8ba07609340ea6">total_num_pixels</a> +         <span class="comment">// Offset to corresponding scatter image</span>
<a name="l00547"></a>00547                     (pixel_coord_x + pixel_coord_z*(detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>)) ),  <span class="comment">// Offset to the corresponding pixel</span>
<a name="l00548"></a>00548                     __float2ull_rn((*energy)*<a class="code" href="MC-GPU__v1_83_8h.html#ac0c03cf5e803404f5dc2a42467395a3a" title="Value to scale the deposited energy in the pixels so that it can be stored as a long long integer ins...">SCALE_eV</a>) );     <span class="comment">// Energy arriving at the pixel, scaled by the factor SCALE_eV and rounded.</span>
<a name="l00549"></a>00549                                                               <span class="comment">// The maximum unsigned long long int value is ~1.8e19:</span>
<a name="l00550"></a>00550         }
<a name="l00551"></a>00551       }
<a name="l00552"></a>00552 <span class="preprocessor">    #else</span>
<a name="l00553"></a>00553 <span class="preprocessor"></span>      <span class="comment">// CPU version (not using CUDA intrinsics: atomicAdd, fast type casting)</span>
<a name="l00554"></a>00554       rotated_position = detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[0]*position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[1]*position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[2]*position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;  <span class="comment">// X coordinate</span>
<a name="l00555"></a>00555       
<a name="l00556"></a>00556       <span class="keywordtype">float</span> pixel_coord_x = floor((rotated_position - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a3c6f9e59e89db06211ccea63287e2826">inv_pixel_size_X</a>);   <span class="comment">// Using float+floor instead of INT to avoid truncation errors for positive and negative values</span>
<a name="l00557"></a>00557       <span class="keywordflow">if</span> ( (pixel_coord_x&gt;-0.1f) &amp;&amp; (pixel_coord_x&lt;(detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>-0.1f)) )    <span class="comment">// Rejecting values negative or bigger than the image size</span>
<a name="l00558"></a>00558       {
<a name="l00559"></a>00559         rotated_position = detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[6]*position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[7]*position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> + detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a05b03d920370a5925ea56932e5953a69">rot_inv</a>[8]*position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;  <span class="comment">// Z coordinate</span>
<a name="l00560"></a>00560         <span class="keywordtype">float</span> pixel_coord_z = floor((rotated_position - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#aa91fc94f674ad7af8eca577f4d8cd22b">inv_pixel_size_Z</a>);
<a name="l00561"></a>00561         <span class="keywordflow">if</span> ( (pixel_coord_z&gt;-0.1f) &amp;&amp; (pixel_coord_z&lt;(detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>-0.1f)) )
<a name="l00562"></a>00562           image[(int)(((<span class="keywordtype">float</span>)*scatter_state)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#ac26d0b66ebf521346d8ba07609340ea6">total_num_pixels</a> + pixel_coord_x + pixel_coord_z*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>  +  0.0001f)]
<a name="l00563"></a>00563              += (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span>)((*energy)*<a class="code" href="MC-GPU__v1_83_8h.html#ac0c03cf5e803404f5dc2a42467395a3a" title="Value to scale the deposited energy in the pixels so that it can be stored as a long long integer ins...">SCALE_eV</a> + 0.5f);   <span class="comment">// Tally the particle energy in the pixel. This instruction is not thread-safe, but it is ok in sequential CPU code.          </span>
<a name="l00564"></a>00564       }
<a name="l00565"></a>00565 <span class="preprocessor">    #endif</span>
<a name="l00566"></a>00566 <span class="preprocessor"></span>  }
<a name="l00567"></a>00567   <span class="keywordflow">else</span>  <span class="comment">// (detector_data_SHARED-&gt;rotation_flag != 1) --&gt;  Initial source direction is (0,1,0): pixel number and distance can be found easily</span>
<a name="l00568"></a>00568   {  
<a name="l00569"></a>00569     <span class="keywordflow">if</span> (direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &lt; 0.0001f)
<a name="l00570"></a>00570       <span class="keywordflow">return</span>;  <span class="comment">// *** Reject particles not moving towards the detector plane at +Y.</span>
<a name="l00571"></a>00571 
<a name="l00572"></a>00572     dist_detector = (detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1f0ce65af62c45925dfc96de3b796474">center</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> - position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>)/(direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>);  <span class="comment">// Distance to the intersection with the detector at +Y.</span>
<a name="l00573"></a>00573   
<a name="l00574"></a>00574       <span class="comment">// !!DeBuG!! IF below (v1.2) not needed when checking the angle</span>
<a name="l00575"></a>00575       <span class="comment">//     if (dist_detector&gt;(2.1f*detector_data_SHARED-&gt;sdd)) return;  </span>
<a name="l00576"></a>00576      
<a name="l00577"></a>00577     
<a name="l00578"></a>00578 <span class="preprocessor">    #ifdef USING_CUDA</span>
<a name="l00579"></a>00579 <span class="preprocessor"></span>    <span class="keywordtype">int</span> pixel_coord_x = __float2int_rd((position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + dist_detector*direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a3c6f9e59e89db06211ccea63287e2826">inv_pixel_size_X</a>);
<a name="l00580"></a>00580     <span class="keywordflow">if</span> ((pixel_coord_x&gt;-1)&amp;&amp;(pixel_coord_x&lt;detector_data_SHARED-&gt;num_pixels.x))
<a name="l00581"></a>00581     {
<a name="l00582"></a>00582       <span class="keywordtype">int</span> pixel_coord_z = __float2int_rd((position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> + dist_detector*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#aa91fc94f674ad7af8eca577f4d8cd22b">inv_pixel_size_Z</a>);
<a name="l00583"></a>00583       <span class="keywordflow">if</span> ((pixel_coord_z&gt;-1)&amp;&amp;(pixel_coord_z&lt;detector_data_SHARED-&gt;num_pixels.y))
<a name="l00584"></a>00584         atomicAdd( ( image +                                                                <span class="comment">// Pointer to beginning of image array</span>
<a name="l00585"></a>00585                      (<span class="keywordtype">int</span>)(*scatter_state) * detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#ac26d0b66ebf521346d8ba07609340ea6">total_num_pixels</a> +         <span class="comment">// Offset to corresponding scatter image</span>
<a name="l00586"></a>00586                      (pixel_coord_x + pixel_coord_z*(detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>)) ),  <span class="comment">// Offset to the corresponding pixel</span>
<a name="l00587"></a>00587                    __float2ull_rn((*energy)*<a class="code" href="MC-GPU__v1_83_8h.html#ac0c03cf5e803404f5dc2a42467395a3a" title="Value to scale the deposited energy in the pixels so that it can be stored as a long long integer ins...">SCALE_eV</a>) );    <span class="comment">// Energy arriving at the pixel, scaled by the factor SCALE_eV and rounded.</span>
<a name="l00588"></a>00588     }
<a name="l00589"></a>00589 <span class="preprocessor">    #else</span>
<a name="l00590"></a>00590 <span class="preprocessor"></span>
<a name="l00591"></a>00591     <span class="comment">// --Calculate the pixel the xray enters, truncating towards minus infinite and making sure the conversion to int is safe:</span>
<a name="l00592"></a>00592     <span class="keywordtype">float</span> pixel_coord_x = floor((position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + dist_detector*direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a3c6f9e59e89db06211ccea63287e2826">inv_pixel_size_X</a>);
<a name="l00593"></a>00593 
<a name="l00594"></a>00594     <span class="keywordflow">if</span> ( (pixel_coord_x&gt;-0.1f) &amp;&amp; (pixel_coord_x&lt;(detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>-0.1f)) )
<a name="l00595"></a>00595     {
<a name="l00596"></a>00596       <span class="keywordtype">float</span> pixel_coord_z = floor((position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> + dist_detector*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> - detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a824e406500b0d8ca59ca96985a599b8b">corner_min_rotated_to_Y</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#aa91fc94f674ad7af8eca577f4d8cd22b">inv_pixel_size_Z</a>);
<a name="l00597"></a>00597       <span class="keywordflow">if</span> ( (pixel_coord_z&gt;-0.1f) &amp;&amp; (pixel_coord_z&lt;(detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>-0.1f)) )
<a name="l00598"></a>00598         image[(int)(((<span class="keywordtype">float</span>)*scatter_state)*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#ac26d0b66ebf521346d8ba07609340ea6">total_num_pixels</a> + pixel_coord_x + pixel_coord_z*detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#a1a30e7168cba7eebe918ab39bb408dce">num_pixels</a>.<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>  +  0.0001f)]
<a name="l00599"></a>00599            += (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span>)((*energy)*<a class="code" href="MC-GPU__v1_83_8h.html#ac0c03cf5e803404f5dc2a42467395a3a" title="Value to scale the deposited energy in the pixels so that it can be stored as a long long integer ins...">SCALE_eV</a> + 0.5f);    <span class="comment">// Truncate the pixel number to INT and round the energy value</span>
<a name="l00600"></a>00600     }
<a name="l00601"></a>00601 <span class="preprocessor">    #endif</span>
<a name="l00602"></a>00602 <span class="preprocessor"></span>  }
<a name="l00603"></a>00603 
<a name="l00604"></a>00604 }
<a name="l00605"></a>00605 
<a name="l00606"></a>00606 
<a name="l00607"></a>00607 <span class="comment"></span>
<a name="l00608"></a>00608 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00609"></a>00609 <span class="comment">//!  Source that creates primary x rays, according to the defined source model.</span>
<a name="l00610"></a>00610 <span class="comment">//!  The particles are automatically moved to the surface of the voxel bounding box,</span>
<a name="l00611"></a>00611 <span class="comment">//!  to start the tracking inside a real material. If the sampled particle do not</span>
<a name="l00612"></a>00612 <span class="comment">//!  enter the voxels, it is init in the focal spot and the main program will check</span>
<a name="l00613"></a>00613 <span class="comment">//!  if it arrives at the detector or not.</span>
<a name="l00614"></a>00614 <span class="comment">//!</span>
<a name="l00615"></a>00615 <span class="comment">//!       @param[in] source_data   Structure describing the source.</span>
<a name="l00616"></a>00616 <span class="comment">//!       @param[in] source_energy_data_CONST   Global variable in constant memory space describing the source energy spectrum.</span>
<a name="l00617"></a>00617 <span class="comment">//!       @param[out] position   Initial particle position (particle transported inside the voxel bbox).</span>
<a name="l00618"></a>00618 <span class="comment">//!       @param[out] direction   Sampled particle direction (cosine vectors).</span>
<a name="l00619"></a>00619 <span class="comment">//!       @param[out] energy   Sampled energy of the new x ray.</span>
<a name="l00620"></a>00620 <span class="comment">//!       @param[in] seed   Current seed of the random number generator, requiered to sample the movement direction.</span>
<a name="l00621"></a>00621 <span class="comment">//!       @param[out] absvox   Set to &lt;0 if primary particle will not cross the voxels, not changed otherwise (&gt;0).</span>
<a name="l00622"></a>00622 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00623"></a>00623 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00624"></a>00624 <span class="preprocessor"></span>__device__
<a name="l00625"></a>00625 <span class="preprocessor">#endif</span>
<a name="l00626"></a><a class="code" href="MC-GPU__v1_83_8h.html#a6807df36f6f2501a9f57ea4fce81b4f2">00626</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a6807df36f6f2501a9f57ea4fce81b4f2" title="Source that creates primary x rays, according to the defined source model.">source</a>(<a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>, <a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>, <span class="keywordtype">float</span>* energy, <a class="code" href="structint2.html">int2</a>* seed, <span class="keywordtype">int</span>* absvox, <span class="keyword">struct</span> <a class="code" href="structsource__struct.html" title="Structure storing the data defining the source model (except for the energy spectrum).">source_struct</a>* source_data_SHARED, <span class="keyword">struct</span> <a class="code" href="structdetector__struct.html" title="Structure storing the data defining the x-ray detector.">detector_struct</a>* detector_data_SHARED)
<a name="l00627"></a>00627 {
<a name="l00628"></a>00628   <span class="comment">// *** Sample the initial x-ray energy following the input energy spectrum using the Walker aliasing algorithm from PENELOPE:</span>
<a name="l00629"></a>00629       <span class="comment">// The following code is equivalent to calling the function &quot;seeki_walker&quot;: int sampled_bin = seeki_walker(source_data_CONST.espc_cutoff, source_data_CONST.espc_alias, ranecu(seed), source_data_CONST.num_bins_espc);      </span>
<a name="l00630"></a>00630   <span class="keywordtype">int</span> sampled_bin;
<a name="l00631"></a>00631   <span class="keywordtype">float</span> RN = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed) * <a class="code" href="MC-GPU__v1_83_8h.html#ad18b13e5db1d8c7c4acdd9af492d6914" title="Global variable to be stored in the GPU constant memory defining the source energy spectrum...">source_energy_data_CONST</a>.<a class="code" href="structsource__energy__struct.html#a559e19adbc83b2aea219281e14f82bb1">num_bins_espc</a>;    <span class="comment">// Find initial interval (array starting at 0):   </span>
<a name="l00632"></a>00632 <span class="preprocessor">  #ifdef USING_CUDA</span>
<a name="l00633"></a>00633 <span class="preprocessor"></span>    <span class="keywordtype">int</span> int_part = __float2int_rd(RN);                          <span class="comment">//   -- Integer part (round down)</span>
<a name="l00634"></a>00634 <span class="preprocessor">  #else</span>
<a name="l00635"></a>00635 <span class="preprocessor"></span>    <span class="keywordtype">int</span> int_part = (int)(RN);
<a name="l00636"></a>00636 <span class="preprocessor">  #endif</span>
<a name="l00637"></a>00637 <span class="preprocessor"></span>  <span class="keywordtype">float</span> fraction_part = RN - ((float)int_part);                 <span class="comment">//   -- Fractional part</span>
<a name="l00638"></a>00638   <span class="keywordflow">if</span> (fraction_part &lt; <a class="code" href="MC-GPU__v1_83_8h.html#ad18b13e5db1d8c7c4acdd9af492d6914" title="Global variable to be stored in the GPU constant memory defining the source energy spectrum...">source_energy_data_CONST</a>.<a class="code" href="structsource__energy__struct.html#a39635503e804b265f76d3a0f7edb2a9f">espc_cutoff</a>[int_part])  <span class="comment">// Check if we are in the aliased part</span>
<a name="l00639"></a>00639     sampled_bin = int_part;                                     <span class="comment">// Below the cutoff: return current value</span>
<a name="l00640"></a>00640   <span class="keywordflow">else</span>
<a name="l00641"></a>00641     sampled_bin = (int)<a class="code" href="MC-GPU__v1_83_8h.html#ad18b13e5db1d8c7c4acdd9af492d6914" title="Global variable to be stored in the GPU constant memory defining the source energy spectrum...">source_energy_data_CONST</a>.<a class="code" href="structsource__energy__struct.html#a8504b5d435cbf4feef2af89c0fd3f829">espc_alias</a>[int_part];  <span class="comment">// Above the cutoff: return alias</span>
<a name="l00642"></a>00642   
<a name="l00643"></a>00643   <span class="comment">// Linear interpolation of the final energy within the sampled energy bin:</span>
<a name="l00644"></a>00644   *energy = <a class="code" href="MC-GPU__v1_83_8h.html#ad18b13e5db1d8c7c4acdd9af492d6914" title="Global variable to be stored in the GPU constant memory defining the source energy spectrum...">source_energy_data_CONST</a>.<a class="code" href="structsource__energy__struct.html#a66db20f7a335eb91023e0d1191730b42">espc</a>[sampled_bin] + <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed) * (<a class="code" href="MC-GPU__v1_83_8h.html#ad18b13e5db1d8c7c4acdd9af492d6914" title="Global variable to be stored in the GPU constant memory defining the source energy spectrum...">source_energy_data_CONST</a>.<a class="code" href="structsource__energy__struct.html#a66db20f7a335eb91023e0d1191730b42">espc</a>[sampled_bin+1] - <a class="code" href="MC-GPU__v1_83_8h.html#ad18b13e5db1d8c7c4acdd9af492d6914" title="Global variable to be stored in the GPU constant memory defining the source energy spectrum...">source_energy_data_CONST</a>.<a class="code" href="structsource__energy__struct.html#a66db20f7a335eb91023e0d1191730b42">espc</a>[sampled_bin]);   
<a name="l00645"></a>00645       
<a name="l00646"></a>00646  
<a name="l00647"></a>00647    <span class="comment">// *** Sample the initial direction:</span>
<a name="l00648"></a>00648    
<a name="l00649"></a>00649   <span class="keywordflow">do</span>   <span class="comment">//  Iterate sampling if the sampled direction is not acceptable to get a square field at the given phi (rejection sampling): force square field for any phi!!</span>
<a name="l00650"></a>00650   {
<a name="l00651"></a>00651     <span class="comment">//     Using the algorithm used in PENMAIN.f, from penelope 2008 (by F. Salvat).</span>
<a name="l00652"></a>00652     direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ad3b7a93633ec87d036453e54e7e7abad">cos_theta_low</a> + <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed)*source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a942289f702f8fc1c6acd411c8ff19efe">D_cos_theta</a>;     <span class="comment">// direction-&gt;z = w = cos(theta_sampled)</span>
<a name="l00653"></a>00653     <span class="keyword">register</span> <span class="keywordtype">float</span> phi_sampled = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a22bb360d294d1a3236e8ba8c4b6d06aa">phi_low</a> + <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed)*source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a01b06b35c2bc440ad0091dddde0cad06">D_phi</a>;
<a name="l00654"></a>00654     <span class="keyword">register</span> <span class="keywordtype">float</span> sin_theta_sampled = sqrtf(1.0f - direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>);
<a name="l00655"></a>00655     <span class="keywordtype">float</span> sinphi_sampled, cosphi_sampled;
<a name="l00656"></a>00656     
<a name="l00657"></a>00657 <span class="preprocessor">    #ifdef USING_CUDA</span>
<a name="l00658"></a>00658 <span class="preprocessor"></span>      sincos(phi_sampled, &amp;sinphi_sampled,&amp;cosphi_sampled);    <span class="comment">// Calculate the SIN and COS at the same time.</span>
<a name="l00659"></a>00659 <span class="preprocessor">    #else</span>
<a name="l00660"></a>00660 <span class="preprocessor"></span>      sinphi_sampled = sin(phi_sampled);   <span class="comment">// Some CPU compilers will be able to use &quot;sincos&quot;, but let&#39;s be safe.</span>
<a name="l00661"></a>00661       cosphi_sampled = cos(phi_sampled);
<a name="l00662"></a>00662 <span class="preprocessor">    #endif       </span>
<a name="l00663"></a>00663 <span class="preprocessor"></span>    
<a name="l00664"></a>00664     direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = sin_theta_sampled * sinphi_sampled;
<a name="l00665"></a>00665     direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = sin_theta_sampled * cosphi_sampled;
<a name="l00666"></a>00666   }
<a name="l00667"></a>00667   <span class="keywordflow">while</span>( fabsf(direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>/(direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>+1.0e-7f)) &gt; source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#ab1bdb59150b1b055e9149461bc090cf4">max_height_at_y1cm</a> );  <span class="comment">// !!DeBuG!! Force square field for any phi by rejection sampling!! Is it necessary to use the &quot;+1.0e-7f&quot; to prevent possible division by zero???</span>
<a name="l00668"></a>00668     
<a name="l00669"></a>00669 
<a name="l00670"></a>00670   <span class="keywordflow">if</span> (detector_data_SHARED-&gt;<a class="code" href="structdetector__struct.html#adbe41a01886ae097dac0fa908551d035">rotation_flag</a> == 1)
<a name="l00671"></a>00671   {
<a name="l00672"></a>00672     <span class="comment">// --Initial beam not pointing to (0,1,0), apply rotation:</span>
<a name="l00673"></a>00673     <span class="keyword">register</span> <span class="keywordtype">float</span> direction_x_tmp = direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l00674"></a>00674     <span class="keyword">register</span> <span class="keywordtype">float</span> direction_y_tmp = direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l00675"></a>00675     direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[0]*direction_x_tmp + source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[1]*direction_y_tmp + source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[2]*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00676"></a>00676     direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[3]*direction_x_tmp + source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[4]*direction_y_tmp + source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[5]*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00677"></a>00677     direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[6]*direction_x_tmp + source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[7]*direction_y_tmp + source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#a0fcd590a93652a70c73fbc8577d36cfb">rot_fan</a>[8]*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00678"></a>00678   }
<a name="l00679"></a>00679 
<a name="l00680"></a>00680   <span class="comment">// Initialize x ray position at focal spot before translation into bbox. Particle stays in focal spot if no interaction found:</span>
<a name="l00681"></a>00681   position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l00682"></a>00682   position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l00683"></a>00683   position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = source_data_SHARED-&gt;<a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00684"></a>00684       
<a name="l00685"></a>00685   <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a047931bac84f3bce50dc13061371e765" title="Functions that moves a particle inside the voxelized geometry bounding box.">move_to_bbox</a>(position, direction, <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#abdd12181cdff77da8b880eedf1d81e33">size_bbox</a>, absvox);  <span class="comment">// Move the particle inside the voxel bounding box.</span>
<a name="l00686"></a>00686 }
<a name="l00687"></a>00687 
<a name="l00688"></a>00688 
<a name="l00689"></a>00689 <span class="comment"></span>
<a name="l00690"></a>00690 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00691"></a>00691 <span class="comment">//!  Functions that moves a particle inside the voxelized geometry bounding box.</span>
<a name="l00692"></a>00692 <span class="comment">//!  An EPSILON distance is added to make sure the particles will be clearly inside the bbox, </span>
<a name="l00693"></a>00693 <span class="comment">//!  not exactly on the surface. </span>
<a name="l00694"></a>00694 <span class="comment">//!</span>
<a name="l00695"></a>00695 <span class="comment">//!  This algorithm makes the following assumtions:</span>
<a name="l00696"></a>00696 <span class="comment">//!     - The back lower vertex of the voxel bounding box is always located at the origin: (x0,y0,z0)=(0,0,0).</span>
<a name="l00697"></a>00697 <span class="comment">//!     - The initial value of &quot;position&quot; corresponds to the focal spot location.</span>
<a name="l00698"></a>00698 <span class="comment">//!     - When a ray is not pointing towards the bbox plane that it should cross according to the sign of the direction,</span>
<a name="l00699"></a>00699 <span class="comment">//!       I assign a distance to the intersection =0 instead of the real negative distance. The wall that will be </span>
<a name="l00700"></a>00700 <span class="comment">//!       crossed to enter the bbox is always the furthest and therefore a 0 distance will never be used except</span>
<a name="l00701"></a>00701 <span class="comment">//!       in the case of a ray starting inside the bbox or outside the bbox and not pointing to any of the 3 planes. </span>
<a name="l00702"></a>00702 <span class="comment">//!       In this situation the ray will be transported a 0 distance, meaning that it will stay at the focal spot.</span>
<a name="l00703"></a>00703 <span class="comment">//!</span>
<a name="l00704"></a>00704 <span class="comment">//!  (Interesting information on ray-box intersection: http://tog.acm.org/resources/GraphicsGems/gems/RayBox.c)</span>
<a name="l00705"></a>00705 <span class="comment">//!</span>
<a name="l00706"></a>00706 <span class="comment">//!       @param[in,out] position Particle position: initially set to the focal spot, returned transported inside the voxel bbox.</span>
<a name="l00707"></a>00707 <span class="comment">//!       @param[out] direction   Sampled particle direction (cosine vectors).</span>
<a name="l00708"></a>00708 <span class="comment">//!       @param[out] intersection_flag   Set to &lt;0 if particle outside bbox and will not cross the voxels, not changed otherwise.</span>
<a name="l00709"></a>00709 <span class="comment">//!       @param[out] size_bbox   Size of the bounding box.</span>
<a name="l00710"></a>00710 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00711"></a>00711 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00712"></a>00712 <span class="preprocessor"></span>__device__
<a name="l00713"></a>00713 <span class="preprocessor">#endif</span>
<a name="l00714"></a><a class="code" href="MC-GPU__v1_83_8h.html#a047931bac84f3bce50dc13061371e765">00714</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a047931bac84f3bce50dc13061371e765" title="Functions that moves a particle inside the voxelized geometry bounding box.">move_to_bbox</a>(<a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>, <a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>, <a class="code" href="structfloat3.html">float3</a> size_bbox, <span class="keywordtype">int</span>* intersection_flag)
<a name="l00715"></a>00715 {
<a name="l00716"></a>00716   <span class="keywordtype">float</span> dist_y, dist_x, dist_z;
<a name="l00717"></a>00717 
<a name="l00718"></a>00718   <span class="comment">// -Distance to the nearest Y plane:</span>
<a name="l00719"></a>00719   <span class="keywordflow">if</span> ((direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>) &gt; <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>)   <span class="comment">// Moving to +Y: check distance to y=0 plane</span>
<a name="l00720"></a>00720   {
<a name="l00721"></a>00721     <span class="comment">// Check Y=0 (bbox wall):</span>
<a name="l00722"></a>00722     <span class="keywordflow">if</span> (position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &gt; 0.0f)  <span class="comment">// The input position must correspond to the focal spot =&gt; position-&gt;y == source_data_CONST.position[*num_p].y</span>
<a name="l00723"></a>00723       dist_y = 0.0f;  <span class="comment">// No intersection with this plane: particle inside or past the box  </span>
<a name="l00724"></a>00724           <span class="comment">// The actual distance would be negative but we set it to 0 bc we will not move the particle if no intersection exist.</span>
<a name="l00725"></a>00725     <span class="keywordflow">else</span>
<a name="l00726"></a>00726       dist_y = <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a> + (-position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>)/(direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>);    <span class="comment">// dist_y &gt; 0 for sure in this case</span>
<a name="l00727"></a>00727   }
<a name="l00728"></a>00728   <span class="keywordflow">else</span> <span class="keywordflow">if</span> ((direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>) &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a0bc5a0759d92c51db359dddd471d6537">NEG_EPS_SOURCE</a>)
<a name="l00729"></a>00729   {
<a name="l00730"></a>00730     <span class="comment">// Check Y=voxel_data_CONST.size_bbox.y:</span>
<a name="l00731"></a>00731     <span class="keywordflow">if</span> (position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &lt; size_bbox.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>)
<a name="l00732"></a>00732       dist_y = 0.0f;  <span class="comment">// No intersection with this plane</span>
<a name="l00733"></a>00733     <span class="keywordflow">else</span>
<a name="l00734"></a>00734       dist_y = <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a> + (size_bbox.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> - position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>)/(direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>);    <span class="comment">// dist_y &gt; 0 for sure in this case</span>
<a name="l00735"></a>00735   }
<a name="l00736"></a>00736   <span class="keywordflow">else</span>   <span class="comment">// (direction-&gt;y)~0</span>
<a name="l00737"></a>00737     dist_y = <a class="code" href="MC-GPU__v1_83_8h.html#a98e3e165debaf418b6b285ff7538d3f3">NEG_INF</a>;   <span class="comment">// Particle moving parallel to the plane: no interaction possible (set impossible negative dist = -INFINITE)</span>
<a name="l00738"></a>00738 
<a name="l00739"></a>00739   <span class="comment">// -Distance to the nearest X plane:</span>
<a name="l00740"></a>00740   <span class="keywordflow">if</span> ((direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>) &gt; <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>)
<a name="l00741"></a>00741   {
<a name="l00742"></a>00742     <span class="comment">// Check X=0:</span>
<a name="l00743"></a>00743     <span class="keywordflow">if</span> (position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> &gt; 0.0f)
<a name="l00744"></a>00744       dist_x = 0.0f;
<a name="l00745"></a>00745     <span class="keywordflow">else</span>  
<a name="l00746"></a>00746       dist_x = <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a> + (-position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>)/(direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>);    <span class="comment">// dist_x &gt; 0 for sure in this case</span>
<a name="l00747"></a>00747   }
<a name="l00748"></a>00748   <span class="keywordflow">else</span> <span class="keywordflow">if</span> ((direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>) &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a0bc5a0759d92c51db359dddd471d6537">NEG_EPS_SOURCE</a>)
<a name="l00749"></a>00749   {
<a name="l00750"></a>00750     <span class="comment">// Check X=voxel_data_CONST.size_bbox.x:</span>
<a name="l00751"></a>00751     <span class="keywordflow">if</span> (position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> &lt; size_bbox.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>)
<a name="l00752"></a>00752       dist_x = 0.0f;
<a name="l00753"></a>00753     <span class="keywordflow">else</span>  
<a name="l00754"></a>00754       dist_x = <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a> + (size_bbox.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> - position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>)/(direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>);    <span class="comment">// dist_x &gt; 0 for sure in this case</span>
<a name="l00755"></a>00755   }
<a name="l00756"></a>00756   <span class="keywordflow">else</span>
<a name="l00757"></a>00757     dist_x = <a class="code" href="MC-GPU__v1_83_8h.html#a98e3e165debaf418b6b285ff7538d3f3">NEG_INF</a>;
<a name="l00758"></a>00758 
<a name="l00759"></a>00759   <span class="comment">// -Distance to the nearest Z plane:</span>
<a name="l00760"></a>00760   <span class="keywordflow">if</span> ((direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>) &gt; <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>)
<a name="l00761"></a>00761   {
<a name="l00762"></a>00762     <span class="comment">// Check Z=0:</span>
<a name="l00763"></a>00763     <span class="keywordflow">if</span> (position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> &gt; 0.0f)
<a name="l00764"></a>00764       dist_z = 0.0f;
<a name="l00765"></a>00765     <span class="keywordflow">else</span>
<a name="l00766"></a>00766       dist_z = <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a> + (-position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>)/(direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>);    <span class="comment">// dist_z &gt; 0 for sure in this case</span>
<a name="l00767"></a>00767   }
<a name="l00768"></a>00768   <span class="keywordflow">else</span> <span class="keywordflow">if</span> ((direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>) &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a0bc5a0759d92c51db359dddd471d6537">NEG_EPS_SOURCE</a>)
<a name="l00769"></a>00769   {
<a name="l00770"></a>00770     <span class="comment">// Check Z=voxel_data_CONST.size_bbox.z:</span>
<a name="l00771"></a>00771     <span class="keywordflow">if</span> (position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> &lt; size_bbox.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>)
<a name="l00772"></a>00772       dist_z = 0.0f;
<a name="l00773"></a>00773     <span class="keywordflow">else</span>
<a name="l00774"></a>00774       dist_z = <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a> + (size_bbox.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> - position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>)/(direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>);    <span class="comment">// dist_z &gt; 0 for sure in this case</span>
<a name="l00775"></a>00775   }
<a name="l00776"></a>00776   <span class="keywordflow">else</span>
<a name="l00777"></a>00777     dist_z = <a class="code" href="MC-GPU__v1_83_8h.html#a98e3e165debaf418b6b285ff7538d3f3">NEG_INF</a>;
<a name="l00778"></a>00778 
<a name="l00779"></a>00779   
<a name="l00780"></a>00780   <span class="comment">// -- Find the longest distance plane, which is the one that has to be crossed to enter the bbox.</span>
<a name="l00781"></a>00781   <span class="comment">//    Storing the maximum distance in variable &quot;dist_z&quot;. Distance will be =0 if no intersection exists or </span>
<a name="l00782"></a>00782   <span class="comment">//    if the x ray is already inside the bbox.</span>
<a name="l00783"></a>00783   <span class="keywordflow">if</span> ( (dist_y&gt;dist_x) &amp;&amp; (dist_y&gt;dist_z) )
<a name="l00784"></a>00784     dist_z = dist_y;      <span class="comment">// dist_z == dist_max </span>
<a name="l00785"></a>00785   <span class="keywordflow">else</span> <span class="keywordflow">if</span> (dist_x&gt;dist_z)
<a name="l00786"></a>00786     dist_z = dist_x;
<a name="l00787"></a>00787 <span class="comment">// else</span>
<a name="l00788"></a>00788 <span class="comment">//   dist_max = dist_z;</span>
<a name="l00789"></a>00789     
<a name="l00790"></a>00790   <span class="comment">// -- Move particle from the focal spot (current location) to the bbox wall surface (slightly inside):</span>
<a name="l00791"></a>00791   position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> += dist_z * direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l00792"></a>00792   position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> += dist_z * direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l00793"></a>00793   position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> += dist_z * direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;      
<a name="l00794"></a>00794   
<a name="l00795"></a>00795   <span class="comment">// Check if the new position is outside the bbox. If true, the particle must be moved back to the focal spot location:</span>
<a name="l00796"></a>00796   <span class="keywordflow">if</span> ( (position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> &lt; 0.0f) || (position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> &gt; size_bbox.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>) || 
<a name="l00797"></a>00797        (position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &lt; 0.0f) || (position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &gt; size_bbox.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>) || 
<a name="l00798"></a>00798        (position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> &lt; 0.0f) || (position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> &gt; size_bbox.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>) )
<a name="l00799"></a>00799   {
<a name="l00800"></a>00800     position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> -= dist_z * direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;  <span class="comment">// Reject new position undoing the previous translation</span>
<a name="l00801"></a>00801     position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> -= dist_z * direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l00802"></a>00802     position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> -= dist_z * direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l00803"></a>00803     (*intersection_flag) = -111;  <span class="comment">// Particle outside the bbox AND not pointing to the bbox: set absvox&lt;0 to skip interaction sampling.</span>
<a name="l00804"></a>00804   }
<a name="l00805"></a>00805 }
<a name="l00806"></a>00806 
<a name="l00807"></a>00807 <span class="comment"></span>
<a name="l00808"></a>00808 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00809"></a>00809 <span class="comment"></span>
<a name="l00810"></a>00810 <span class="comment"></span>
<a name="l00811"></a>00811 <span class="comment">//!  Upper limit of the number of random values sampled in a single track.</span>
<a name="l00812"></a><a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3ab75c5bf72afd4a9b0255f61bd95704">00812</a> <span class="comment"></span><span class="preprocessor">#define  LEAP_DISTANCE     256</span>
<a name="l00813"></a>00813 <span class="preprocessor"></span><span class="comment">//!  Multipliers and moduli for the two MLCG in RANECU.</span>
<a name="l00814"></a><a class="code" href="MC-GPU__kernel__v1_83_8cu.html#afd7d5320ad79ad99c8d230f3a730141d">00814</a> <span class="comment"></span><span class="preprocessor">#define  a1_RANECU       40014</span>
<a name="l00815"></a><a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aa7d2cfcd3830e91736ca79f4bfdfb5d2">00815</a> <span class="preprocessor"></span><span class="preprocessor">#define  m1_RANECU  2147483563</span>
<a name="l00816"></a><a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a527a42aca15fac6a1cab1116a72eb8ca">00816</a> <span class="preprocessor"></span><span class="preprocessor">#define  a2_RANECU       40692</span>
<a name="l00817"></a><a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a9872bd70bdcc27f56bc025f291b7d6c2">00817</a> <span class="preprocessor"></span><span class="preprocessor">#define  m2_RANECU  2147483399</span>
<a name="l00818"></a>00818 <span class="preprocessor"></span><span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00819"></a>00819 <span class="comment"></span><span class="comment">//! Initialize the pseudo-random number generator (PRNG) RANECU to a position</span>
<a name="l00820"></a>00820 <span class="comment"></span><span class="comment">//! far away from the previous history (leap frog technique).</span>
<a name="l00821"></a>00821 <span class="comment"></span><span class="comment">//!</span>
<a name="l00822"></a>00822 <span class="comment"></span><span class="comment">//! Each calculated seed initiates a consecutive and disjoint sequence of</span>
<a name="l00823"></a>00823 <span class="comment"></span><span class="comment">//! pseudo-random numbers with length LEAP_DISTANCE, that can be used to</span>
<a name="l00824"></a>00824 <span class="comment"></span><span class="comment">//! in a parallel simulation (Sequence Splitting parallelization method).</span>
<a name="l00825"></a>00825 <span class="comment"></span><span class="comment">//! The basic equation behind the algorithm is:</span>
<a name="l00826"></a>00826 <span class="comment"></span><span class="comment">//!    S(i+j) = (a**j * S(i)) MOD m = [(a**j MOD m)*S(i)] MOD m  ,</span>
<a name="l00827"></a>00827 <span class="comment"></span><span class="comment">//! which is described in:</span>
<a name="l00828"></a>00828 <span class="comment"></span><span class="comment">//!   P L&#39;Ecuyer, Commun. ACM 31 (1988) p.742</span>
<a name="l00829"></a>00829 <span class="comment"></span><span class="comment">//!</span>
<a name="l00830"></a>00830 <span class="comment"></span><span class="comment">//! This function has been adapted from &quot;seedsMLCG.f&quot;, see:</span>
<a name="l00831"></a>00831 <span class="comment"></span><span class="comment">//!   A Badal and J Sempau, Computer Physics Communications 175 (2006) p. 440-450</span>
<a name="l00832"></a>00832 <span class="comment"></span><span class="comment">//!</span>
<a name="l00833"></a>00833 <span class="comment"></span><span class="comment">//!       @param[in] history   Particle bach number.</span>
<a name="l00834"></a>00834 <span class="comment"></span><span class="comment">//!       @param[in] seed_input   Initial PRNG seed input (used to initiate both MLCGs in RANECU).</span>
<a name="l00835"></a>00835 <span class="comment"></span><span class="comment">//!       @param[out] seed   Initial PRNG seeds for the present history.</span>
<a name="l00836"></a>00836 <span class="comment"></span><span class="comment">//!</span>
<a name="l00837"></a>00837 <span class="comment"></span><span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00838"></a>00838 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00839"></a>00839 <span class="preprocessor"></span>__device__
<a name="l00840"></a>00840 <span class="preprocessor">#endif</span>
<a name="l00841"></a><a class="code" href="MC-GPU__v1_83_8h.html#ae7bd0b93b4d5eedb15822fa0ef622b48">00841</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae7bd0b93b4d5eedb15822fa0ef622b48" title="Initialize the pseudo-random number generator (PRNG) RANECU to a position far away from the previous ...">init_PRNG</a>(<span class="keywordtype">int</span> history_batch, <span class="keywordtype">int</span> histories_per_thread, <span class="keywordtype">int</span> seed_input, <a class="code" href="structint2.html">int2</a>* seed)
<a name="l00842"></a>00842 {
<a name="l00843"></a>00843   <span class="comment">// -- Move the RANECU generator to a unique position for the current batch of histories:</span>
<a name="l00844"></a>00844   <span class="comment">//    I have to use an &quot;unsigned long long int&quot; value to represent all the simulated histories in all previous batches</span>
<a name="l00845"></a>00845   <span class="comment">//    The maximum unsigned long long int value is ~1.8e19: if history &gt;1.8e16 and LEAP_DISTANCE==1000, &#39;leap&#39; will overflow.</span>
<a name="l00846"></a>00846   <span class="comment">// **** 1st MLCG:</span>
<a name="l00847"></a>00847   <span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> <span class="keywordtype">int</span> leap = ((<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)(history_batch+1))*(histories_per_thread*<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3ab75c5bf72afd4a9b0255f61bd95704" title="Upper limit of the number of random values sampled in a single track.">LEAP_DISTANCE</a>);
<a name="l00848"></a>00848   <span class="keywordtype">int</span> y = 1;
<a name="l00849"></a>00849   <span class="keywordtype">int</span> z = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#afd7d5320ad79ad99c8d230f3a730141d" title="Multipliers and moduli for the two MLCG in RANECU.">a1_RANECU</a>;
<a name="l00850"></a>00850   <span class="comment">// -- Calculate the modulo power &#39;(a^leap)MOD(m)&#39; using a divide-and-conquer algorithm adapted to modulo arithmetic</span>
<a name="l00851"></a>00851   <span class="keywordflow">for</span>(;;)
<a name="l00852"></a>00852   {
<a name="l00853"></a>00853     <span class="comment">// (A2) Halve n, and store the integer part and the residue</span>
<a name="l00854"></a>00854     <span class="keywordflow">if</span> (0!=(leap&amp;01))  <span class="comment">// (bit-wise operation for MOD(leap,2), or leap%2 ==&gt; proceed if leap is an odd number)  Equivalent: t=(short)(leap%2);</span>
<a name="l00855"></a>00855     {
<a name="l00856"></a>00856       leap &gt;&gt;= 1;     <span class="comment">// Halve n moving the bits 1 position right. Equivalent to:  leap=(leap/2);  </span>
<a name="l00857"></a>00857       y = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aa7d2cfcd3830e91736ca79f4bfdfb5d2">m1_RANECU</a>,z,y);      <span class="comment">// (A3) Multiply y by z:  y = [z*y] MOD m</span>
<a name="l00858"></a>00858       <span class="keywordflow">if</span> (0==leap) <span class="keywordflow">break</span>;         <span class="comment">// (A4) leap==0? ==&gt; finish</span>
<a name="l00859"></a>00859     }
<a name="l00860"></a>00860     <span class="keywordflow">else</span>           <span class="comment">// (leap is even)</span>
<a name="l00861"></a>00861     {
<a name="l00862"></a>00862       leap&gt;&gt;= 1;     <span class="comment">// Halve leap moving the bits 1 position right. Equivalent to:  leap=(leap/2);</span>
<a name="l00863"></a>00863     }
<a name="l00864"></a>00864     z = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aa7d2cfcd3830e91736ca79f4bfdfb5d2">m1_RANECU</a>,z,z);        <span class="comment">// (A5) Square z:  z = [z*z] MOD m</span>
<a name="l00865"></a>00865   }
<a name="l00866"></a>00866   <span class="comment">// AjMODm1 = y;                 // Exponentiation finished:  AjMODm = expMOD = y = a^j</span>
<a name="l00867"></a>00867 
<a name="l00868"></a>00868   <span class="comment">// -- Compute and display the seeds S(i+j), from the present seed S(i), using the previously calculated value of (a^j)MOD(m):</span>
<a name="l00869"></a>00869   <span class="comment">//         S(i+j) = [(a**j MOD m)*S(i)] MOD m</span>
<a name="l00870"></a>00870   <span class="comment">//         S_i = abMODm(m,S_i,AjMODm)</span>
<a name="l00871"></a>00871   seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#aa7d2cfcd3830e91736ca79f4bfdfb5d2">m1_RANECU</a>, seed_input, y);     <span class="comment">// Using the input seed as the starting seed</span>
<a name="l00872"></a>00872 
<a name="l00873"></a>00873   <span class="comment">// **** 2nd MLCG (repeating the previous calculation for the 2nd MLCG parameters):</span>
<a name="l00874"></a>00874   leap = ((<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)(history_batch+1))*(histories_per_thread*LEAP_DISTANCE);
<a name="l00875"></a>00875   y = 1;
<a name="l00876"></a>00876   z = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a527a42aca15fac6a1cab1116a72eb8ca">a2_RANECU</a>;
<a name="l00877"></a>00877   <span class="keywordflow">for</span>(;;)
<a name="l00878"></a>00878   {
<a name="l00879"></a>00879     <span class="comment">// (A2) Halve n, and store the integer part and the residue</span>
<a name="l00880"></a>00880     <span class="keywordflow">if</span> (0!=(leap&amp;01))  <span class="comment">// (bit-wise operation for MOD(leap,2), or leap%2 ==&gt; proceed if leap is an odd number)  Equivalent: t=(short)(leap%2);</span>
<a name="l00881"></a>00881     {
<a name="l00882"></a>00882       leap &gt;&gt;= 1;     <span class="comment">// Halve n moving the bits 1 position right. Equivalent to:  leap=(leap/2);</span>
<a name="l00883"></a>00883       y = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a9872bd70bdcc27f56bc025f291b7d6c2">m2_RANECU</a>,z,y);      <span class="comment">// (A3) Multiply y by z:  y = [z*y] MOD m</span>
<a name="l00884"></a>00884       <span class="keywordflow">if</span> (0==leap) <span class="keywordflow">break</span>;         <span class="comment">// (A4) leap==0? ==&gt; finish</span>
<a name="l00885"></a>00885     }
<a name="l00886"></a>00886     <span class="keywordflow">else</span>           <span class="comment">// (leap is even)</span>
<a name="l00887"></a>00887     {
<a name="l00888"></a>00888       leap&gt;&gt;= 1;     <span class="comment">// Halve leap moving the bits 1 position right. Equivalent to:  leap=(leap/2);</span>
<a name="l00889"></a>00889     }
<a name="l00890"></a>00890     z = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a9872bd70bdcc27f56bc025f291b7d6c2">m2_RANECU</a>,z,z);        <span class="comment">// (A5) Square z:  z = [z*z] MOD m</span>
<a name="l00891"></a>00891   }
<a name="l00892"></a>00892   <span class="comment">// AjMODm2 = y;</span>
<a name="l00893"></a>00893   seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a9872bd70bdcc27f56bc025f291b7d6c2">m2_RANECU</a>, seed_input, y);     <span class="comment">// Using the input seed as the starting seed</span>
<a name="l00894"></a>00894 }
<a name="l00895"></a>00895 
<a name="l00896"></a>00896 
<a name="l00897"></a>00897 <span class="comment"></span>
<a name="l00898"></a>00898 <span class="comment">/////////////////////////////////////////////////////////////////////</span>
<a name="l00899"></a>00899 <span class="comment">//!  Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding</span>
<a name="l00900"></a>00900 <span class="comment">//!  the possible overflow, using the Russian Peasant approach</span>
<a name="l00901"></a>00901 <span class="comment">//!  modulo m and the approximate factoring method, as described</span>
<a name="l00902"></a>00902 <span class="comment">//!  in:  L&#39;Ecuyer and Cote, ACM Trans. Math. Soft. 17 (1991).</span>
<a name="l00903"></a>00903 <span class="comment">//!</span>
<a name="l00904"></a>00904 <span class="comment">//!  This function has been adapted from &quot;seedsMLCG.f&quot;, see: </span>
<a name="l00905"></a>00905 <span class="comment">//!  Badal and Sempau, Computer Physics Communications 175 (2006)</span>
<a name="l00906"></a>00906 <span class="comment">//!</span>
<a name="l00907"></a>00907 <span class="comment">//!       @param[in] m,a,s  MLCG parameters</span>
<a name="l00908"></a>00908 <span class="comment">//!       @return   (a1*a2) MOD m   </span>
<a name="l00909"></a>00909 <span class="comment"></span><span class="comment">//</span>
<a name="l00910"></a>00910 <span class="comment">//    Input:          0 &lt; a1 &lt; m                                  </span>
<a name="l00911"></a>00911 <span class="comment">//                    0 &lt; a2 &lt; m                                  </span>
<a name="l00912"></a>00912 <span class="comment">//</span>
<a name="l00913"></a>00913 <span class="comment">//    Return value:  (a1*a2) MOD m                                </span>
<a name="l00914"></a>00914 <span class="comment">//</span><span class="comment"></span>
<a name="l00915"></a>00915 <span class="comment">/////////////////////////////////////////////////////////////////////</span>
<a name="l00916"></a>00916 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00917"></a>00917 <span class="preprocessor"></span>__device__ __host__    <span class="comment">// Function will be callable from host and also from device</span>
<a name="l00918"></a>00918 <span class="preprocessor">#endif</span>
<a name="l00919"></a><a class="code" href="MC-GPU__v1_83_8h.html#ae9061f4314474bf4415200a5cb115397">00919</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">int</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#ae9061f4314474bf4415200a5cb115397" title="Calculate &quot;(a1*a2) MOD m&quot; with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L&#39;Ecuyer and Cote, ACM Trans.">abMODm</a>(<span class="keywordtype">int</span> m, <span class="keywordtype">int</span> a, <span class="keywordtype">int</span> s)
<a name="l00920"></a>00920 {
<a name="l00921"></a>00921   <span class="comment">// CAUTION: the input parameters are modified in the function but should not be returned to the calling function! (pass by value!)</span>
<a name="l00922"></a>00922   <span class="keywordtype">int</span> q, k;
<a name="l00923"></a>00923   <span class="keywordtype">int</span> p = -m;            <span class="comment">// p is always negative to avoid overflow when adding</span>
<a name="l00924"></a>00924 
<a name="l00925"></a>00925   <span class="comment">// ** Apply the Russian peasant method until &quot;a =&lt; 32768&quot;:</span>
<a name="l00926"></a>00926   <span class="keywordflow">while</span> (a&gt;32768)        <span class="comment">// We assume &#39;32&#39; bit integers (4 bytes): 2^((&#39;32&#39;-2)/2) = 32768</span>
<a name="l00927"></a>00927   {
<a name="l00928"></a>00928     <span class="keywordflow">if</span> (0!=(a&amp;1))        <span class="comment">// Store &#39;s&#39; when &#39;a&#39; is odd     Equivalent code:   if (1==(a%2))</span>
<a name="l00929"></a>00929     {
<a name="l00930"></a>00930       p += s;
<a name="l00931"></a>00931       <span class="keywordflow">if</span> (p&gt;0) p -= m;
<a name="l00932"></a>00932     }
<a name="l00933"></a>00933     a &gt;&gt;= 1;             <span class="comment">// Half a (move bits 1 position right)   Equivalent code: a = a/2;</span>
<a name="l00934"></a>00934     s = (s-m) + s;       <span class="comment">// Double s (MOD m)</span>
<a name="l00935"></a>00935     <span class="keywordflow">if</span> (s&lt;0) s += m;     <span class="comment">// (s is always positive)</span>
<a name="l00936"></a>00936   }
<a name="l00937"></a>00937 
<a name="l00938"></a>00938   <span class="comment">// ** Employ the approximate factoring method (a is small enough to avoid overflow):</span>
<a name="l00939"></a>00939   q = (int) m / a;
<a name="l00940"></a>00940   k = (int) s / q;
<a name="l00941"></a>00941   s = a*(s-k*q)-k*(m-q*a);
<a name="l00942"></a>00942   <span class="keywordflow">while</span> (s&lt;0)
<a name="l00943"></a>00943     s += m;
<a name="l00944"></a>00944 
<a name="l00945"></a>00945   <span class="comment">// ** Compute the final result:</span>
<a name="l00946"></a>00946   p += s;
<a name="l00947"></a>00947   <span class="keywordflow">if</span> (p&lt;0) p += m;
<a name="l00948"></a>00948 
<a name="l00949"></a>00949   <span class="keywordflow">return</span> p;
<a name="l00950"></a>00950 }
<a name="l00951"></a>00951 
<a name="l00952"></a>00952 
<a name="l00953"></a>00953 <span class="comment"></span>
<a name="l00954"></a>00954 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00955"></a>00955 <span class="comment">//! Pseudo-random number generator (PRNG) RANECU returning a float value</span>
<a name="l00956"></a>00956 <span class="comment">//! (single precision version).</span>
<a name="l00957"></a>00957 <span class="comment">//!</span>
<a name="l00958"></a>00958 <span class="comment">//!       @param[in,out] seed   PRNG seed (seed kept in the calling function and updated here).</span>
<a name="l00959"></a>00959 <span class="comment">//!       @return   PRN double value in the open interval (0,1)</span>
<a name="l00960"></a>00960 <span class="comment">//!</span>
<a name="l00961"></a>00961 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00962"></a>00962 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00963"></a>00963 <span class="preprocessor"></span>__device__ 
<a name="l00964"></a>00964 <span class="preprocessor">#endif</span>
<a name="l00965"></a><a class="code" href="MC-GPU__v1_83_8h.html#a3b57f2faa5765ceba5d352641ca72a9f">00965</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">float</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(<a class="code" href="structint2.html">int2</a>* seed)
<a name="l00966"></a>00966 {
<a name="l00967"></a>00967   <span class="keywordtype">int</span> i1 = (int)(seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>/53668);
<a name="l00968"></a>00968   seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> = 40014*(seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>-i1*53668)-i1*12211;
<a name="l00969"></a>00969 
<a name="l00970"></a>00970   <span class="keywordtype">int</span> i2 = (int)(seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>/52774);
<a name="l00971"></a>00971   seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> = 40692*(seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>-i2*52774)-i2*3791;
<a name="l00972"></a>00972 
<a name="l00973"></a>00973   <span class="keywordflow">if</span> (seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> &lt; 0) seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> += 2147483563;
<a name="l00974"></a>00974   <span class="keywordflow">if</span> (seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> &lt; 0) seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> += 2147483399;
<a name="l00975"></a>00975 
<a name="l00976"></a>00976   i2 = seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>-seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>;
<a name="l00977"></a>00977   <span class="keywordflow">if</span> (i2 &lt; 1) i2 += 2147483562;
<a name="l00978"></a>00978 
<a name="l00979"></a>00979 
<a name="l00980"></a>00980 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00981"></a>00981 <span class="preprocessor"></span>  <span class="keywordflow">return</span> (__int2float_rn(i2)*4.65661305739e-10f);        <span class="comment">// 4.65661305739e-10 == 1/2147483563</span>
<a name="l00982"></a>00982 <span class="preprocessor">#else</span>
<a name="l00983"></a>00983 <span class="preprocessor"></span>  <span class="keywordflow">return</span> ((<span class="keywordtype">float</span>)(i2)*4.65661305739e-10f);          
<a name="l00984"></a>00984 <span class="preprocessor">#endif</span>
<a name="l00985"></a>00985 <span class="preprocessor"></span>
<a name="l00986"></a>00986 }
<a name="l00987"></a>00987 
<a name="l00988"></a>00988 <span class="comment"></span>
<a name="l00989"></a>00989 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00990"></a>00990 <span class="comment">//! Pseudo-random number generator (PRNG) RANECU returning a double value.</span>
<a name="l00991"></a>00991 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l00992"></a>00992 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l00993"></a>00993 <span class="preprocessor"></span>__device__ 
<a name="l00994"></a>00994 <span class="preprocessor">#endif</span>
<a name="l00995"></a><a class="code" href="MC-GPU__v1_83_8h.html#a74c4894502584f8fb113e04b759a6c55">00995</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">double</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(<a class="code" href="structint2.html">int2</a>* seed)
<a name="l00996"></a>00996 {
<a name="l00997"></a>00997   <span class="keywordtype">int</span> i1 = (int)(seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>/53668);
<a name="l00998"></a>00998   seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> = 40014*(seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>-i1*53668)-i1*12211;
<a name="l00999"></a>00999 
<a name="l01000"></a>01000   <span class="keywordtype">int</span> i2 = (int)(seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>/52774);
<a name="l01001"></a>01001   seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> = 40692*(seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>-i2*52774)-i2*3791;
<a name="l01002"></a>01002 
<a name="l01003"></a>01003   <span class="keywordflow">if</span> (seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> &lt; 0) seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a> += 2147483563;
<a name="l01004"></a>01004   <span class="keywordflow">if</span> (seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> &lt; 0) seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a> += 2147483399;
<a name="l01005"></a>01005 
<a name="l01006"></a>01006   i2 = seed-&gt;<a class="code" href="structint2.html#a2b4908d6f7ae421393c0ec0233e90551">x</a>-seed-&gt;<a class="code" href="structint2.html#a3d7ff3bd7734f8c2f38613a4da5e6c83">y</a>;
<a name="l01007"></a>01007   <span class="keywordflow">if</span> (i2 &lt; 1) i2 += 2147483562;
<a name="l01008"></a>01008 
<a name="l01009"></a>01009 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01010"></a>01010 <span class="preprocessor"></span>  <span class="keywordflow">return</span> (__int2double_rn(i2)*4.6566130573917692e-10);
<a name="l01011"></a>01011 <span class="preprocessor">#else</span>
<a name="l01012"></a>01012 <span class="preprocessor"></span>  <span class="keywordflow">return</span> ((<span class="keywordtype">double</span>)(i2)*4.6566130573917692e-10);
<a name="l01013"></a>01013 <span class="preprocessor">#endif</span>
<a name="l01014"></a>01014 <span class="preprocessor"></span>
<a name="l01015"></a>01015 }
<a name="l01016"></a>01016 
<a name="l01017"></a>01017 
<a name="l01018"></a>01018 <span class="comment"></span>
<a name="l01019"></a>01019 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l01020"></a>01020 <span class="comment">//! Find the voxel that contains the current position.</span>
<a name="l01021"></a>01021 <span class="comment">//! Report the voxel absolute index and the x,y,z indices.</span>
<a name="l01022"></a>01022 <span class="comment">//! The structure containing the voxel number and size is read from CONSTANT memory.</span>
<a name="l01023"></a>01023 <span class="comment">//!</span>
<a name="l01024"></a>01024 <span class="comment">//!       @param[in] position   Particle position</span>
<a name="l01025"></a>01025 <span class="comment">//!       @param[out] voxel_coord   Pointer to three integer values (short3*) that will store the x,y and z voxel indices.</span>
<a name="l01026"></a>01026 <span class="comment">//!       @return   Returns &quot;absvox&quot;, the voxel number where the particle is</span>
<a name="l01027"></a>01027 <span class="comment">//!                 located (negative if position outside the voxel bbox).</span>
<a name="l01028"></a>01028 <span class="comment">//!</span>
<a name="l01029"></a>01029 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l01030"></a>01030 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01031"></a>01031 <span class="preprocessor"></span>__device__
<a name="l01032"></a>01032 <span class="preprocessor">#endif</span>
<a name="l01033"></a><a class="code" href="MC-GPU__v1_83_8h.html#a562dba755f1abdb20b1088d3121f4326">01033</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">int</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a562dba755f1abdb20b1088d3121f4326" title="Find the voxel that contains the current position.">locate_voxel</a>(<a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#af0d7a559e56da9df6b01dda2ac95a787">position</a>, <a class="code" href="structshort3.html">short3</a>* voxel_coord)
<a name="l01034"></a>01034 {
<a name="l01035"></a>01035 
<a name="l01036"></a>01036   <span class="keywordflow">if</span> ( (position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>) || (position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> &gt; (<a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#abdd12181cdff77da8b880eedf1d81e33">size_bbox</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> - <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>)) ||
<a name="l01037"></a>01037        (position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>) || (position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> &gt; (<a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#abdd12181cdff77da8b880eedf1d81e33">size_bbox</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> - <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>)) ||
<a name="l01038"></a>01038        (position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> &lt; <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>) || (position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> &gt; (<a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#abdd12181cdff77da8b880eedf1d81e33">size_bbox</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> - <a class="code" href="MC-GPU__v1_83_8h.html#a4b98b60caf9dba36fca9fe41c3238404" title="Offset value to make sure the particles are completely inside, or outside, the voxel bounding box...">EPS_SOURCE</a>)) )
<a name="l01039"></a>01039   {
<a name="l01040"></a>01040     <span class="comment">// -- Particle escaped the voxelized geometry (using EPS_SOURCE to avoid numerical precision errors):      </span>
<a name="l01041"></a>01041      <span class="keywordflow">return</span> -1;
<a name="l01042"></a>01042   }
<a name="l01043"></a>01043  
<a name="l01044"></a>01044   <span class="comment">// -- Particle inside the voxelized geometry, find current voxel:</span>
<a name="l01045"></a>01045   <span class="comment">//    The truncation from float to integer could give troubles for negative coordinates but this will never happen thanks to the IF at the begining of this function.</span>
<a name="l01046"></a>01046   <span class="comment">//    (no need to use the CUDA function to convert float to integer rounding down (towards minus infinite): __float2int_rd)</span>
<a name="l01047"></a>01047   
<a name="l01048"></a>01048   <span class="keyword">register</span> <span class="keywordtype">int</span> voxel_coord_x, voxel_coord_y, voxel_coord_z;
<a name="l01049"></a>01049 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01050"></a>01050 <span class="preprocessor"></span>  voxel_coord_x = __float2int_rd(position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> * <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#ad1468b63b2b99455904cc43e9058d5c1">inv_voxel_size</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>);  
<a name="l01051"></a>01051   voxel_coord_y = __float2int_rd(position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> * <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#ad1468b63b2b99455904cc43e9058d5c1">inv_voxel_size</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>);
<a name="l01052"></a>01052   voxel_coord_z = __float2int_rd(position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> * <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#ad1468b63b2b99455904cc43e9058d5c1">inv_voxel_size</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>);
<a name="l01053"></a>01053 <span class="preprocessor">#else</span>
<a name="l01054"></a>01054 <span class="preprocessor"></span>  voxel_coord_x = (int)(position-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> * <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#ad1468b63b2b99455904cc43e9058d5c1">inv_voxel_size</a>.<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>);     
<a name="l01055"></a>01055   voxel_coord_y = (int)(position-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> * <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#ad1468b63b2b99455904cc43e9058d5c1">inv_voxel_size</a>.<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>);
<a name="l01056"></a>01056   voxel_coord_z = (int)(position-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> * <a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#ad1468b63b2b99455904cc43e9058d5c1">inv_voxel_size</a>.<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>);
<a name="l01057"></a>01057 <span class="preprocessor">#endif</span>
<a name="l01058"></a>01058 <span class="preprocessor"></span>
<a name="l01059"></a>01059   <span class="comment">// Output the voxel coordinates as short int (2 bytes) instead of int (4 bytes) to save registers; avoid type castings in the calculation of the return value.</span>
<a name="l01060"></a>01060   voxel_coord-&gt;<a class="code" href="structshort3.html#a5d5d88cacca4de50a2334be0033a126a">x</a> = (<span class="keywordtype">short</span> int) voxel_coord_x;
<a name="l01061"></a>01061   voxel_coord-&gt;<a class="code" href="structshort3.html#a6056b3f7d064291aaf9121575ec83866">y</a> = (<span class="keywordtype">short</span> int) voxel_coord_y;
<a name="l01062"></a>01062   voxel_coord-&gt;<a class="code" href="structshort3.html#a803658b690ee5726c60a2b99b99b9143">z</a> = (<span class="keywordtype">short</span> int) voxel_coord_z;
<a name="l01063"></a>01063   
<a name="l01064"></a>01064   <span class="keywordflow">return</span> (voxel_coord_x + voxel_coord_y*(<a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#adeed35f7e44a7cfb2891001a1a68f868">num_voxels</a>.<a class="code" href="structint3.html#a0a4ad50a155a35fa938ce6f16930affa">x</a>) + voxel_coord_z*(<a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#adeed35f7e44a7cfb2891001a1a68f868">num_voxels</a>.<a class="code" href="structint3.html#a0a4ad50a155a35fa938ce6f16930affa">x</a>)*(<a class="code" href="MC-GPU__v1_83_8h.html#a129396d7723d11e10c673c1171773cdc" title="Global variable to be stored in the GPU constant memory defining the size of the voxel phantom...">voxel_data_CONST</a>.<a class="code" href="structvoxel__struct.html#adeed35f7e44a7cfb2891001a1a68f868">num_voxels</a>.<a class="code" href="structint3.html#a5d95e23491677d61019f0354b16adca9">y</a>));  
<a name="l01065"></a>01065 }
<a name="l01066"></a>01066 
<a name="l01067"></a>01067 
<a name="l01068"></a>01068 <span class="comment"></span>
<a name="l01069"></a>01069 <span class="comment">//////////////////////////////////////////////////////////////////////</span>
<a name="l01070"></a>01070 <span class="comment">//!   Rotates a vector; the rotation is specified by giving</span>
<a name="l01071"></a>01071 <span class="comment">//!   the polar and azimuthal angles in the &quot;self-frame&quot;, as</span>
<a name="l01072"></a>01072 <span class="comment">//!   determined by the vector to be rotated.</span>
<a name="l01073"></a>01073 <span class="comment">//!   This function is a literal translation from Fortran to C of</span>
<a name="l01074"></a>01074 <span class="comment">//!   PENELOPE (v. 2006) subroutine &quot;DIRECT&quot;.</span>
<a name="l01075"></a>01075 <span class="comment">//!</span>
<a name="l01076"></a>01076 <span class="comment">//!    @param[in,out]  (u,v,w)  input vector (=d) in the lab. frame; returns the rotated vector components in the lab. frame</span>
<a name="l01077"></a>01077 <span class="comment">//!    @param[in]  costh  cos(theta), angle between d before and after turn</span>
<a name="l01078"></a>01078 <span class="comment">//!    @param[in]  phi  azimuthal angle (rad) turned by d in its self-frame</span>
<a name="l01079"></a>01079 <span class="comment"></span><span class="comment">//</span>
<a name="l01080"></a>01080 <span class="comment">//    Output:</span>
<a name="l01081"></a>01081 <span class="comment">//      (u,v,w) -&gt; rotated vector components in the lab. frame</span>
<a name="l01082"></a>01082 <span class="comment">//</span>
<a name="l01083"></a>01083 <span class="comment">//    Comments:</span>
<a name="l01084"></a>01084 <span class="comment">//      -&gt; (u,v,w) should have norm=1 on input; if not, it is</span>
<a name="l01085"></a>01085 <span class="comment">//         renormalized on output, provided norm&gt;0.</span>
<a name="l01086"></a>01086 <span class="comment">//      -&gt; The algorithm is based on considering the turned vector</span>
<a name="l01087"></a>01087 <span class="comment">//         d&#39; expressed in the self-frame S&#39;,</span>
<a name="l01088"></a>01088 <span class="comment">//           d&#39; = (sin(th)cos(ph), sin(th)sin(ph), cos(th))</span>
<a name="l01089"></a>01089 <span class="comment">//         and then apply a change of frame from S&#39; to the lab</span>
<a name="l01090"></a>01090 <span class="comment">//         frame. S&#39; is defined as having its z&#39; axis coincident</span>
<a name="l01091"></a>01091 <span class="comment">//         with d, its y&#39; axis perpendicular to z and z&#39; and its</span>
<a name="l01092"></a>01092 <span class="comment">//         x&#39; axis equal to y&#39;*z&#39;. The matrix of the change is then</span>
<a name="l01093"></a>01093 <span class="comment">//                   / uv/rho    -v/rho    u \</span>
<a name="l01094"></a>01094 <span class="comment">//          S -&gt;lab: | vw/rho     u/rho    v |  , rho=(u^2+v^2)^0.5</span>
<a name="l01095"></a>01095 <span class="comment">//                   \ -rho       0        w /</span>
<a name="l01096"></a>01096 <span class="comment">//      -&gt; When rho=0 (w=1 or -1) z and z&#39; are parallel and the y&#39;</span>
<a name="l01097"></a>01097 <span class="comment">//         axis cannot be defined in this way. Instead y&#39; is set to</span>
<a name="l01098"></a>01098 <span class="comment">//         y and therefore either x&#39;=x (if w=1) or x&#39;=-x (w=-1)</span><span class="comment"></span>
<a name="l01099"></a>01099 <span class="comment">//////////////////////////////////////////////////////////////////////</span>
<a name="l01100"></a>01100 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01101"></a>01101 <span class="preprocessor"></span>__device__
<a name="l01102"></a>01102 <span class="preprocessor">#endif</span>
<a name="l01103"></a><a class="code" href="MC-GPU__v1_83_8h.html#acc5186c685730d287039c737445bbe43">01103</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a94ad9ed794eb8fd696e7e357ae391765" title="Rotates a vector; the rotation is specified by giving the polar and azimuthal angles in the &quot;self-fra...">rotate_double</a>(<a class="code" href="structfloat3.html">float3</a>* <a class="code" href="structsource__struct.html#ae1a18388dff97d101b36aeaba9bd363a">direction</a>, <span class="keywordtype">double</span> costh, <span class="keywordtype">double</span> phi)   <span class="comment">// !!DeBuG!! The direction vector is single precision but the rotation is performed in doule precision for increased accuracy.</span>
<a name="l01104"></a>01104 {
<a name="l01105"></a>01105   <span class="keywordtype">double</span> DXY, NORM, cosphi, sinphi, SDT;
<a name="l01106"></a>01106   DXY = direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>*direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>*direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l01107"></a>01107   
<a name="l01108"></a>01108 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01109"></a>01109 <span class="preprocessor"></span>  sincos(phi, &amp;sinphi,&amp;cosphi);   <span class="comment">// Calculate the SIN and COS at the same time.</span>
<a name="l01110"></a>01110 <span class="preprocessor">#else</span>
<a name="l01111"></a>01111 <span class="preprocessor"></span>  sinphi = sin(phi);   <span class="comment">// Some CPU compilers will be able to use &quot;sincos&quot;, but let&#39;s be safe.</span>
<a name="l01112"></a>01112   cosphi = cos(phi);
<a name="l01113"></a>01113 <span class="preprocessor">#endif   </span>
<a name="l01114"></a>01114 <span class="preprocessor"></span>
<a name="l01115"></a>01115   <span class="comment">// ****  Ensure normalisation</span>
<a name="l01116"></a>01116   NORM = DXY + direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;     <span class="comment">// !!DeBuG!! Check if it is really necessary to renormalize in a real simulation!!</span>
<a name="l01117"></a>01117   <span class="keywordflow">if</span> (fabs(NORM-1.0)&gt;1.0e-14)
<a name="l01118"></a>01118   {
<a name="l01119"></a>01119     NORM = 1.0/sqrt(NORM);
<a name="l01120"></a>01120     direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = NORM*direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l01121"></a>01121     direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = NORM*direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l01122"></a>01122     direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = NORM*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>;
<a name="l01123"></a>01123     DXY = direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>*direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> + direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>*direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>;
<a name="l01124"></a>01124   }
<a name="l01125"></a>01125   <span class="keywordflow">if</span> (DXY&gt;1.0e-28)
<a name="l01126"></a>01126   {
<a name="l01127"></a>01127     SDT = sqrt((1.0-costh*costh)/DXY);
<a name="l01128"></a>01128     <span class="keywordtype">float</span> direction_x_in = direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>;
<a name="l01129"></a>01129     direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a>*costh + SDT*(direction_x_in*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>*cosphi-direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>*sinphi);
<a name="l01130"></a>01130     direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>*costh+SDT*(direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a>*direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>*cosphi+direction_x_in*sinphi);
<a name="l01131"></a>01131     direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>*costh-DXY*SDT*cosphi;
<a name="l01132"></a>01132   }
<a name="l01133"></a>01133   <span class="keywordflow">else</span>
<a name="l01134"></a>01134   {
<a name="l01135"></a>01135     SDT = sqrt(1.0-costh*costh);
<a name="l01136"></a>01136     direction-&gt;<a class="code" href="structfloat3.html#aa6147d421a81889971f8c66aa92abf0d">y</a> = SDT*sinphi;
<a name="l01137"></a>01137     <span class="keywordflow">if</span> (direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a>&gt;0.0)
<a name="l01138"></a>01138     {
<a name="l01139"></a>01139       direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> = SDT*cosphi;
<a name="l01140"></a>01140       direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> = costh;
<a name="l01141"></a>01141     }
<a name="l01142"></a>01142     <span class="keywordflow">else</span>
<a name="l01143"></a>01143     {
<a name="l01144"></a>01144       direction-&gt;<a class="code" href="structfloat3.html#af621f02abb1c788738fe61ea9807ff9c">x</a> =-SDT*cosphi;
<a name="l01145"></a>01145       direction-&gt;<a class="code" href="structfloat3.html#a772dffd42d89f350c5a1b766c4703245">z</a> =-costh;
<a name="l01146"></a>01146     }
<a name="l01147"></a>01147   }
<a name="l01148"></a>01148 }
<a name="l01149"></a>01149 
<a name="l01150"></a>01150 <span class="comment"></span>
<a name="l01151"></a>01151 <span class="comment">//////////////////////////////////////////////////////////////////////</span>
<a name="l01152"></a>01152 <span class="comment"></span>
<a name="l01153"></a>01153 
<a name="l01154"></a>01154 <span class="comment">//  ***********************************************************************</span>
<a name="l01155"></a>01155 <span class="comment">//  *   Translation of PENELOPE&#39;s &quot;SUBROUTINE GRAa&quot; from FORTRAN77 to C   *</span>
<a name="l01156"></a>01156 <span class="comment">//  ***********************************************************************</span><span class="comment"></span>
<a name="l01157"></a>01157 <span class="comment">//!  Sample a Rayleigh interaction using the sampling algorithm</span>
<a name="l01158"></a>01158 <span class="comment">//!  used in PENELOPE 2006.</span>
<a name="l01159"></a>01159 <span class="comment">//!</span>
<a name="l01160"></a>01160 <span class="comment">//!       @param[in] energy   Particle energy (not modified with Rayleigh)</span>
<a name="l01161"></a>01161 <span class="comment">//!       @param[out] costh_Rayleigh   Cosine of the angular deflection</span>
<a name="l01162"></a>01162 <span class="comment">//!       @param[in] material  Current voxel material</span>
<a name="l01163"></a>01163 <span class="comment"></span><span class="comment">//</span>
<a name="l01164"></a>01164 <span class="comment">//  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC</span>
<a name="l01165"></a>01165 <span class="comment">//  C  PENELOPE/PENGEOM (version 2006)                                     C</span>
<a name="l01166"></a>01166 <span class="comment">//  C    Copyright (c) 2001-2006                                           C</span>
<a name="l01167"></a>01167 <span class="comment">//  C    Universitat de Barcelona                                          C</span>
<a name="l01168"></a>01168 <span class="comment">//  C  Permission to use, copy, modify, distribute and sell this software  C</span>
<a name="l01169"></a>01169 <span class="comment">//  C  and its documentation for any purpose is hereby granted without     C</span>
<a name="l01170"></a>01170 <span class="comment">//  C  fee, provided that the above copyright notice appears in all        C</span>
<a name="l01171"></a>01171 <span class="comment">//  C  copies and that both that copyright notice and this permission      C</span>
<a name="l01172"></a>01172 <span class="comment">//  C  notice appear in all supporting documentation. The Universitat de   C</span>
<a name="l01173"></a>01173 <span class="comment">//  C  Barcelona makes no representations about the suitability of this    C</span>
<a name="l01174"></a>01174 <span class="comment">//  C  software for any purpose. It is provided &quot;as is&quot; without express    C</span>
<a name="l01175"></a>01175 <span class="comment">//  C  or implied warranty.                                                C</span>
<a name="l01176"></a>01176 <span class="comment">//  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC</span><span class="comment"></span>
<a name="l01177"></a>01177 <span class="comment">//////////////////////////////////////////////////////////////////////</span>
<a name="l01178"></a>01178 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01179"></a>01179 <span class="preprocessor"></span>__device__
<a name="l01180"></a>01180 <span class="preprocessor">#endif</span>
<a name="l01181"></a><a class="code" href="MC-GPU__v1_83_8h.html#a5cafd4ebb3d155746b16ca7d231ecad1">01181</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a5cafd4ebb3d155746b16ca7d231ecad1" title="Sample a Rayleigh interaction using the sampling algorithm used in PENELOPE 2006.">GRAa</a>(<span class="keywordtype">float</span> *energy, <span class="keywordtype">double</span> *costh_Rayleigh, <span class="keywordtype">int</span> *mat, <span class="keywordtype">float</span> *pmax_current, <a class="code" href="structint2.html">int2</a> *seed, <span class="keyword">struct</span> <a class="code" href="structrayleigh__struct.html" title="Structure storing the data of the Rayleigh interaction sampling model (equivalent to PENELOPE&#39;s commo...">rayleigh_struct</a>* cgra)
<a name="l01182"></a>01182 {
<a name="l01183"></a>01183 <span class="comment">/*  ****  Energy grid and interpolation constants for the current energy. */</span>
<a name="l01184"></a>01184     <span class="keywordtype">double</span>  xmax = ((double)*energy) * 8.065535669099010e-5;       <span class="comment">// 8.065535669099010e-5 == 2.0*20.6074/510998.918</span>
<a name="l01185"></a>01185     <span class="keywordtype">double</span> x2max = <a class="code" href="MC-GPU__v1_83_8h.html#a9fa6cc963655f85859604988fa68a9e2">min_value</a>( (xmax*xmax) , ((<span class="keywordtype">double</span>)cgra-&gt;<a class="code" href="structrayleigh__struct.html#a3d91ed8dc782c166106093ae028f0dd8">xco</a>[(*mat+1)*<a class="code" href="MC-GPU__v1_83_8h.html#a50e5fe2bed6d3ba8c550beb81e5e40e5">NP_RAYLEIGH</a> - 1]) );   <span class="comment">// Get the last tabulated value of xco for this mat</span>
<a name="l01186"></a>01186     
<a name="l01187"></a>01187     <span class="keywordflow">if</span> (xmax &lt; 0.01)
<a name="l01188"></a>01188     {
<a name="l01189"></a>01189        <span class="keywordflow">do</span>
<a name="l01190"></a>01190        {
<a name="l01191"></a>01191           *costh_Rayleigh = 1.0 - <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(seed) * 2.0;
<a name="l01192"></a>01192        }
<a name="l01193"></a>01193        <span class="keywordflow">while</span> ( <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(seed) &gt; (((*costh_Rayleigh)*(*costh_Rayleigh)+1.0)*0.5) );
<a name="l01194"></a>01194        <span class="keywordflow">return</span>;
<a name="l01195"></a>01195     }
<a name="l01196"></a>01196 
<a name="l01197"></a>01197     <span class="keywordflow">for</span>(;;)    <span class="comment">// (Loop will iterate everytime the sampled value is rejected or above maximum)</span>
<a name="l01198"></a>01198     {
<a name="l01199"></a>01199       <span class="keywordtype">double</span> ru = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(seed) * (double)(*pmax_current);    <span class="comment">// Pmax for the current energy is entered as a parameter</span>
<a name="l01200"></a>01200  
<a name="l01201"></a>01201 <span class="comment">/*  ****  Selection of the interval  (binary search within pre-calculated limits). */</span>
<a name="l01202"></a>01202       <span class="keywordtype">int</span> itn = (int)(ru * (<a class="code" href="MC-GPU__v1_83_8h.html#a50e5fe2bed6d3ba8c550beb81e5e40e5">NP_RAYLEIGH</a>-1));     <span class="comment">// &#39;itn&#39; will never reach the last interval &#39;NP_RAYLEIGH-1&#39;, but this is how RITA is implemented in PENELOPE</span>
<a name="l01203"></a>01203       <span class="keywordtype">int</span> i__ = (int)cgra-&gt;<a class="code" href="structrayleigh__struct.html#a81e58bece94cd34133d55960cee9fe73">itlco</a>[itn + (*mat)*<a class="code" href="MC-GPU__v1_83_8h.html#a50e5fe2bed6d3ba8c550beb81e5e40e5">NP_RAYLEIGH</a>];
<a name="l01204"></a>01204       <span class="keywordtype">int</span> j   = (int)cgra-&gt;<a class="code" href="structrayleigh__struct.html#af8eb91b5cbbc3bdf6d84e3aa9525946a">ituco</a>[itn + (*mat)*<a class="code" href="MC-GPU__v1_83_8h.html#a50e5fe2bed6d3ba8c550beb81e5e40e5">NP_RAYLEIGH</a>];
<a name="l01205"></a>01205       
<a name="l01206"></a>01206       <span class="keywordflow">if</span> ((j - i__) &gt; 1)
<a name="l01207"></a>01207       {
<a name="l01208"></a>01208         <span class="keywordflow">do</span>
<a name="l01209"></a>01209         {
<a name="l01210"></a>01210           <span class="keyword">register</span> <span class="keywordtype">int</span> k = (i__ + j)&gt;&gt;1;     <span class="comment">// &gt;&gt;1 == /2 </span>
<a name="l01211"></a>01211           <span class="keywordflow">if</span> (ru &gt; cgra-&gt;<a class="code" href="structrayleigh__struct.html#a85ba59a95287c8d2ec3895b9ea6f8f43">pco</a>[k -1 + (*mat)*<a class="code" href="MC-GPU__v1_83_8h.html#a50e5fe2bed6d3ba8c550beb81e5e40e5">NP_RAYLEIGH</a>])
<a name="l01212"></a>01212             i__ = k;
<a name="l01213"></a>01213           <span class="keywordflow">else</span>
<a name="l01214"></a>01214             j = k;
<a name="l01215"></a>01215         }
<a name="l01216"></a>01216         <span class="keywordflow">while</span> ((j - i__) &gt; 1);
<a name="l01217"></a>01217       }
<a name="l01218"></a>01218        
<a name="l01219"></a>01219 <span class="comment">/*  ****  Sampling from the rational inverse cumulative distribution. */</span>
<a name="l01220"></a>01220       <span class="keywordtype">int</span> index = i__ - 1 + (*mat)*<a class="code" href="MC-GPU__v1_83_8h.html#a50e5fe2bed6d3ba8c550beb81e5e40e5">NP_RAYLEIGH</a>;
<a name="l01221"></a>01221 
<a name="l01222"></a>01222       <span class="keywordtype">double</span> rr = ru - cgra-&gt;<a class="code" href="structrayleigh__struct.html#a85ba59a95287c8d2ec3895b9ea6f8f43">pco</a>[index];
<a name="l01223"></a>01223       <span class="keywordtype">double</span> xx;
<a name="l01224"></a>01224       <span class="keywordflow">if</span> (rr &gt; 1e-16)
<a name="l01225"></a>01225       {      
<a name="l01226"></a>01226         <span class="keywordtype">double</span> d__ = (double)(cgra-&gt;<a class="code" href="structrayleigh__struct.html#a85ba59a95287c8d2ec3895b9ea6f8f43">pco</a>[index+1] - cgra-&gt;<a class="code" href="structrayleigh__struct.html#a85ba59a95287c8d2ec3895b9ea6f8f43">pco</a>[index]);
<a name="l01227"></a>01227         <span class="keywordtype">float</span> aco_index = cgra-&gt;<a class="code" href="structrayleigh__struct.html#a2ea1b62a4e4cc5900c060fed2d95fec5">aco</a>[index], bco_index = cgra-&gt;<a class="code" href="structrayleigh__struct.html#af107a7fbdd989e049ba67b9fc2e614c3">bco</a>[index], xco_index = cgra-&gt;<a class="code" href="structrayleigh__struct.html#a3d91ed8dc782c166106093ae028f0dd8">xco</a>[index];   <span class="comment">// Avoid multiple accesses to the same global variable</span>
<a name="l01228"></a>01228 
<a name="l01229"></a>01229         xx = (double)xco_index + (<span class="keywordtype">double</span>)(aco_index + 1.0f + bco_index)* d__* rr / (d__*d__ + (aco_index*d__ + bco_index*rr) * rr) * (<span class="keywordtype">double</span>)(cgra-&gt;<a class="code" href="structrayleigh__struct.html#a3d91ed8dc782c166106093ae028f0dd8">xco</a>[index+1] - xco_index);
<a name="l01230"></a>01230         
<a name="l01231"></a>01231       }
<a name="l01232"></a>01232       <span class="keywordflow">else</span>
<a name="l01233"></a>01233       {
<a name="l01234"></a>01234         xx = cgra-&gt;<a class="code" href="structrayleigh__struct.html#a3d91ed8dc782c166106093ae028f0dd8">xco</a>[index];
<a name="l01235"></a>01235       }
<a name="l01236"></a>01236       
<a name="l01237"></a>01237       <span class="keywordflow">if</span> (xx &lt; x2max)
<a name="l01238"></a>01238       {
<a name="l01239"></a>01239         <span class="comment">// Sampled value below maximum possible value:</span>
<a name="l01240"></a>01240         *costh_Rayleigh = 1.0 - 2.0 * xx / x2max;   <span class="comment">// !!DeBuG!! costh_Rayleigh in double precision, but not all intermediate steps are!?</span>
<a name="l01241"></a>01241         <span class="comment">/*  ****  Rejection: */</span>    
<a name="l01242"></a>01242         <span class="keywordflow">if</span> (<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a74c4894502584f8fb113e04b759a6c55" title="Pseudo-random number generator (PRNG) RANECU returning a double value.">ranecu_double</a>(seed) &lt; (((*costh_Rayleigh)*(*costh_Rayleigh) + 1.0)*0.5))
<a name="l01243"></a>01243           <span class="keywordflow">break</span>;   <span class="comment">// Sample value not rejected! break loop and return.</span>
<a name="l01244"></a>01244       }
<a name="l01245"></a>01245     }
<a name="l01246"></a>01246 } <span class="comment">/* graa */</span>
<a name="l01247"></a>01247 
<a name="l01248"></a>01248 
<a name="l01249"></a>01249 <span class="comment"></span>
<a name="l01250"></a>01250 <span class="comment">//////////////////////////////////////////////////////////////////////////</span>
<a name="l01251"></a>01251 <span class="comment"></span>
<a name="l01252"></a>01252 
<a name="l01253"></a>01253 <span class="comment">//  ***********************************************************************</span>
<a name="l01254"></a>01254 <span class="comment">//  *   Translation of PENELOPE&#39;s &quot;SUBROUTINE GCOa&quot;  from FORTRAN77 to C  *</span>
<a name="l01255"></a>01255 <span class="comment">//  ********************************************************************* *</span><span class="comment"></span>
<a name="l01256"></a>01256 <span class="comment">//!  Random sampling of incoherent (Compton) scattering of photons, using </span>
<a name="l01257"></a>01257 <span class="comment">//!  the sampling algorithm from PENELOPE 2006:</span>
<a name="l01258"></a>01258 <span class="comment">//!    Relativistic impulse approximation with analytical one-electron Compton profiles</span>
<a name="l01259"></a>01259 <span class="comment"></span>
<a name="l01260"></a>01260 <span class="comment">// !!DeBuG!!  In penelope, Doppler broadening is not used for E greater than 5 MeV.</span>
<a name="l01261"></a>01261 <span class="comment">//            We don&#39;t use it in GPU to reduce the lines of code and prevent using COMMON/compos/ZT(M)</span>
<a name="l01262"></a>01262 <span class="comment"></span>
<a name="l01263"></a>01263 <span class="comment">//!       @param[in,out] energy   incident and final photon energy (eV)</span>
<a name="l01264"></a>01264 <span class="comment">//!       @param[out] costh_Compton   cosine of the polar scattering angle</span>
<a name="l01265"></a>01265 <span class="comment">//!       @param[in] material   Current voxel material</span>
<a name="l01266"></a>01266 <span class="comment">//!       @param[in] seed   RANECU PRNG seed</span>
<a name="l01267"></a>01267 <span class="comment"></span><span class="comment">//</span>
<a name="l01268"></a>01268 <span class="comment">//  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC</span>
<a name="l01269"></a>01269 <span class="comment">//  C  PENELOPE/PENGEOM (version 2006)                                     C</span>
<a name="l01270"></a>01270 <span class="comment">//  C    Copyright (c) 2001-2006                                           C</span>
<a name="l01271"></a>01271 <span class="comment">//  C    Universitat de Barcelona                                          C</span>
<a name="l01272"></a>01272 <span class="comment">//  C  Permission to use, copy, modify, distribute and sell this software  C</span>
<a name="l01273"></a>01273 <span class="comment">//  C  and its documentation for any purpose is hereby granted without     C</span>
<a name="l01274"></a>01274 <span class="comment">//  C  fee, provided that the above copyright notice appears in all        C</span>
<a name="l01275"></a>01275 <span class="comment">//  C  copies and that both that copyright notice and this permission      C</span>
<a name="l01276"></a>01276 <span class="comment">//  C  notice appear in all supporting documentation. The Universitat de   C</span>
<a name="l01277"></a>01277 <span class="comment">//  C  Barcelona makes no representations about the suitability of this    C</span>
<a name="l01278"></a>01278 <span class="comment">//  C  software for any purpose. It is provided &quot;as is&quot; without express    C</span>
<a name="l01279"></a>01279 <span class="comment">//  C  or implied warranty.                                                C</span>
<a name="l01280"></a>01280 <span class="comment">//  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC</span>
<a name="l01281"></a>01281 <span class="comment">//</span>
<a name="l01282"></a>01282 <span class="comment">//  ************************************************************************</span>
<a name="l01283"></a>01283 
<a name="l01284"></a>01284 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01285"></a>01285 <span class="preprocessor"></span>__device__
<a name="l01286"></a>01286 <span class="preprocessor">#endif</span>
<a name="l01287"></a><a class="code" href="MC-GPU__v1_83_8h.html#a6f212862e7fba1553ec9495f05443b2c">01287</a> <span class="preprocessor"></span><span class="keyword">inline</span> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a6f212862e7fba1553ec9495f05443b2c" title="Random sampling of incoherent (Compton) scattering of photons, using the sampling algorithm from PENE...">GCOa</a>(<span class="keywordtype">float</span> *energy, <span class="keywordtype">double</span> *costh_Compton, <span class="keywordtype">int</span> *mat, <a class="code" href="structint2.html">int2</a> *seed, <span class="keyword">struct</span> <a class="code" href="structcompton__struct.html" title="Structure storing the data of the Compton interaction sampling model (equivalent to PENELOPE&#39;s common...">compton_struct</a>* cgco_SHARED)
<a name="l01288"></a>01288 {
<a name="l01289"></a>01289     <span class="keywordtype">float</span> s, a1, s0, af, ek, ek2, ek3, tau, pzomc, taumin;
<a name="l01290"></a>01290     <span class="keywordtype">float</span> rn[<a class="code" href="MC-GPU__v1_83_8h.html#a7fc9275784608be88cd40ac55f4de283">MAX_SHELLS</a>];
<a name="l01291"></a>01291     <span class="keywordtype">double</span> cdt1;
<a name="l01292"></a>01292 
<a name="l01293"></a>01293      <span class="comment">// Some variables used in PENELOPE have been eliminated to save register: float aux, taum2, fpzmax, a, a2, ek1 ,rni, xqc, fpz, pac[MAX_SHELLS];</span>
<a name="l01294"></a>01294 
<a name="l01295"></a>01295     <span class="keywordtype">int</span> i__;
<a name="l01296"></a>01296     <span class="keywordtype">int</span> my_noscco = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#aef81eb007c90d13552ea867912f0d2c3">noscco</a>[*mat];    <span class="comment">// Store the number of oscillators for the input material in a local variable</span>
<a name="l01297"></a>01297     
<a name="l01298"></a>01298 <span class="preprocessor">#ifndef USING_CUDA</span>
<a name="l01299"></a>01299 <span class="preprocessor"></span>    <span class="keyword">static</span> <span class="keywordtype">int</span> warning_flag_1 = -1, warning_flag_2 = -1, warning_flag_3 = -1;    <span class="comment">// Write warnings for the CPU code, but only once.  !!DeBuG!!</span>
<a name="l01300"></a>01300 <span class="preprocessor">#endif</span>
<a name="l01301"></a>01301 <span class="preprocessor"></span>
<a name="l01302"></a>01302     ek = *energy * 1.956951306108245e-6f;    <span class="comment">// (1.956951306108245e-6 == 1.0/510998.918)</span>
<a name="l01303"></a>01303     ek2 = ek * 2.f + 1.f;
<a name="l01304"></a>01304     ek3 = ek * ek;
<a name="l01305"></a>01305     <span class="comment">// ek1 = ek3 - ek2 - 1.;</span>
<a name="l01306"></a>01306     taumin = 1.f / ek2;
<a name="l01307"></a>01307     <span class="comment">// taum2 = taumin * taumin;</span>
<a name="l01308"></a>01308     a1 = logf(ek2);
<a name="l01309"></a>01309     <span class="comment">// a2 = a1 + ek * 2. * (ek + 1.) * taum2;    // a2 was used only once, code moved below</span>
<a name="l01310"></a>01310 
<a name="l01311"></a>01311 
<a name="l01312"></a>01312 <span class="comment">/*  ****  Incoherent scattering function for theta=PI. */</span>
<a name="l01313"></a>01313 
<a name="l01314"></a>01314     s0 = 0.0f;
<a name="l01315"></a>01315     <span class="keywordflow">for</span> (i__ = 0; i__ &lt; my_noscco; i__++)
<a name="l01316"></a>01316     {
<a name="l01317"></a>01317        <span class="keyword">register</span> <span class="keywordtype">float</span> temp = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a32dcebc4af75e0620abed9c32ffc2bdd">uico</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>];
<a name="l01318"></a>01318        <span class="keywordflow">if</span> (temp &lt; *energy)
<a name="l01319"></a>01319        {
<a name="l01320"></a>01320          <span class="keyword">register</span> <span class="keywordtype">float</span> aux = *energy * (*energy - temp) * 2.f;
<a name="l01321"></a>01321 <span class="preprocessor">         #ifdef USING_CUDA</span>
<a name="l01322"></a>01322 <span class="preprocessor"></span>           pzomc = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * (aux - temp * 510998.918f) * rsqrtf(aux + aux + temp * temp) * 1.956951306108245e-6f;
<a name="l01323"></a>01323              <span class="comment">// 1.956951306108245e-6 = 1.0/510998.918f   // Version using the reciprocal of sqrt in CUDA: faster and more accurate!!</span>
<a name="l01324"></a>01324 <span class="preprocessor">         #else</span>
<a name="l01325"></a>01325 <span class="preprocessor"></span>           pzomc = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * (aux - temp * 510998.918f) / (sqrtf(aux + aux + temp * temp) * 510998.918f);
<a name="l01326"></a>01326 <span class="preprocessor">         #endif</span>
<a name="l01327"></a>01327 <span class="preprocessor"></span>         <span class="keywordflow">if</span> (pzomc &gt; 0.0f)
<a name="l01328"></a>01328            temp = (0.707106781186545f+pzomc*1.4142135623731f) * (0.707106781186545f+pzomc*1.4142135623731f);
<a name="l01329"></a>01329          <span class="keywordflow">else</span>
<a name="l01330"></a>01330            temp = (0.707106781186545f-pzomc*1.4142135623731f) * (0.707106781186545f-pzomc*1.4142135623731f);
<a name="l01331"></a>01331 
<a name="l01332"></a>01332          temp = 0.5f * expf(0.5f - temp);    <span class="comment">// Calculate EXP outside the IF to avoid branching</span>
<a name="l01333"></a>01333 
<a name="l01334"></a>01334          <span class="keywordflow">if</span> (pzomc &gt; 0.0f)
<a name="l01335"></a>01335             temp = 1.0f - temp;
<a name="l01336"></a>01336                                 
<a name="l01337"></a>01337          s0 += cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a181131326c64f04da9eb7f199604371e">fco</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * temp;
<a name="l01338"></a>01338        }
<a name="l01339"></a>01339     }
<a name="l01340"></a>01340             
<a name="l01341"></a>01341 <span class="comment">/*  ****  Sampling tau. */</span>
<a name="l01342"></a>01342     <span class="keywordflow">do</span>
<a name="l01343"></a>01343     {
<a name="l01344"></a>01344       <span class="keywordflow">if</span> (<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed)*<span class="comment">/*a2=*/</span>(a1+2.*ek*(ek+1.f)*taumin*taumin) &lt; a1)
<a name="l01345"></a>01345       { 
<a name="l01346"></a>01346         tau = powf(taumin, <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed));    <span class="comment">// !!DeBuG!!  &quot;powf()&quot; has a big error (7 ULP), the double version has only 2!! </span>
<a name="l01347"></a>01347       }
<a name="l01348"></a>01348       <span class="keywordflow">else</span>
<a name="l01349"></a>01349       {
<a name="l01350"></a>01350         tau = sqrtf(1.f + <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed) * (taumin * taumin - 1.f));
<a name="l01351"></a>01351       }
<a name="l01352"></a>01352 
<a name="l01353"></a>01353       cdt1 = (double)(1.f-tau) / (((double)tau)*((double)*energy)*1.956951306108245e-6);    <span class="comment">// !!DeBuG!! The sampled COS will be double precision, but TAU is not!!!</span>
<a name="l01354"></a>01354 
<a name="l01355"></a>01355       <span class="keywordflow">if</span> (cdt1 &gt; 2.0) cdt1 = 1.99999999;   <span class="comment">// !!DeBuG!! Make sure that precision error in POW, SQRT never gives cdt1&gt;2 ==&gt; costh_Compton&lt;-1</span>
<a name="l01356"></a>01356       
<a name="l01357"></a>01357   <span class="comment">/*  ****  Incoherent scattering function. */</span>
<a name="l01358"></a>01358       s = 0.0f;
<a name="l01359"></a>01359       <span class="keywordflow">for</span> (i__ = 0; i__ &lt; my_noscco; i__++)
<a name="l01360"></a>01360       {
<a name="l01361"></a>01361         <span class="keyword">register</span> <span class="keywordtype">float</span> temp = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a32dcebc4af75e0620abed9c32ffc2bdd">uico</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>];
<a name="l01362"></a>01362         <span class="keywordflow">if</span> (temp &lt; *energy)
<a name="l01363"></a>01363         {
<a name="l01364"></a>01364           <span class="keyword">register</span> <span class="keywordtype">float</span> aux = (*energy) * (*energy - temp) * ((<span class="keywordtype">float</span>)cdt1);
<a name="l01365"></a>01365 
<a name="l01366"></a>01366           <span class="keywordflow">if</span> ((aux&gt;1.0e-12f)||(temp&gt;1.0e-12f))  <span class="comment">// !!DeBuG!! Make sure the SQRT argument is never &lt;0, and that we never get 0/0 -&gt; NaN when aux=temp=0 !!</span>
<a name="l01367"></a>01367           {
<a name="l01368"></a>01368 <span class="preprocessor">         #ifdef USING_CUDA</span>
<a name="l01369"></a>01369 <span class="preprocessor"></span>           pzomc = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * (aux - temp * 510998.918f) * rsqrtf(aux + aux + temp * temp) * 1.956951306108245e-6f;
<a name="l01370"></a>01370              <span class="comment">// 1.956951306108245e-6 = 1.0/510998.918f   //  Version using the reciprocal of sqrt in CUDA: faster and more accurate!!</span>
<a name="l01371"></a>01371 <span class="preprocessor">         #else</span>
<a name="l01372"></a>01372 <span class="preprocessor"></span>           pzomc = cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * (aux - temp * 510998.918f) / (sqrtf(aux + aux + temp * temp) * 510998.918f);
<a name="l01373"></a>01373 <span class="preprocessor">         #endif</span>
<a name="l01374"></a>01374 <span class="preprocessor"></span>
<a name="l01375"></a>01375           }
<a name="l01376"></a>01376           <span class="keywordflow">else</span>
<a name="l01377"></a>01377           {
<a name="l01378"></a>01378             pzomc = 0.002f;    <span class="comment">// !!DeBuG!! Using a rough approximation to a sample value of pzomc found using pure double precision: NOT RIGUROUS! But this code is expected to be used very seldom, only in extreme cases.</span>
<a name="l01379"></a>01379 <span class="preprocessor">            #ifndef USING_CUDA</span>
<a name="l01380"></a>01380 <span class="preprocessor"></span>            <span class="keywordflow">if</span> (warning_flag_1&lt;0)
<a name="l01381"></a>01381             {
<a name="l01382"></a>01382                warning_flag_1 = +1;  <span class="comment">// Disable warning, do not show again</span>
<a name="l01383"></a>01383                printf(<span class="stringliteral">&quot;          [... Small numerical precision error detected computing \&quot;pzomc\&quot; in GCOa (this warning will not be repeated).]\n               i__=%d, aux=%.14f, temp=%.14f, pzomc(forced)=%.14f, uico=%.14f, energy=%.7f, cgco_SHARED-&gt;fj0=%.14f, mat=%d, cdt1=%.14lf\n&quot;</span>, (<span class="keywordtype">int</span>)i__, aux, temp, pzomc, cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a32dcebc4af75e0620abed9c32ffc2bdd">uico</a>[*mat+i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>], *energy, cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat+i__*MAX_MATERIALS], (<span class="keywordtype">int</span>)*mat, cdt1);   <span class="comment">// !!DeBuG!!</span>
<a name="l01384"></a>01384             }
<a name="l01385"></a>01385 <span class="preprocessor">            #endif                    </span>
<a name="l01386"></a>01386 <span class="preprocessor"></span>          }
<a name="l01387"></a>01387           
<a name="l01388"></a>01388           temp = pzomc * 1.4142135623731f;
<a name="l01389"></a>01389           <span class="keywordflow">if</span> (pzomc &gt; 0.0f)
<a name="l01390"></a>01390             temp = 0.5f - (temp + 0.70710678118654502f) * (temp + 0.70710678118654502f);   <span class="comment">// Calculate exponential argument</span>
<a name="l01391"></a>01391           <span class="keywordflow">else</span>
<a name="l01392"></a>01392             temp = 0.5f - (0.70710678118654502f - temp) * (0.70710678118654502f - temp);
<a name="l01393"></a>01393 
<a name="l01394"></a>01394           temp = 0.5f * expf(temp);      <span class="comment">// All threads will calculate the expf together</span>
<a name="l01395"></a>01395           
<a name="l01396"></a>01396           <span class="keywordflow">if</span> (pzomc &gt; 0.0f)
<a name="l01397"></a>01397             temp = 1.0f - temp;
<a name="l01398"></a>01398 
<a name="l01399"></a>01399           s += cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a181131326c64f04da9eb7f199604371e">fco</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * temp;
<a name="l01400"></a>01400           rn[i__] = temp;
<a name="l01401"></a>01401         }        
<a name="l01402"></a>01402       }
<a name="l01403"></a>01403     } <span class="keywordflow">while</span>( (<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed)*s0) &gt; (s*(1.0f+tau*(<span class="comment">/*ek1=*/</span>(ek3 - ek2 - 1.0f)+tau*(ek2+tau*ek3)))/(ek3*tau*(tau*tau+1.0f))) );  <span class="comment">//  ****  Rejection function</span>
<a name="l01404"></a>01404 
<a name="l01405"></a>01405     *costh_Compton = 1.0 - cdt1;
<a name="l01406"></a>01406         
<a name="l01407"></a>01407 <span class="comment">/*  ****  Target electron shell. */</span>
<a name="l01408"></a>01408     <span class="keywordflow">for</span> (;;)
<a name="l01409"></a>01409     {
<a name="l01410"></a>01410       <span class="keyword">register</span> <span class="keywordtype">float</span> temp = s*<a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed);
<a name="l01411"></a>01411       <span class="keywordtype">float</span> pac = 0.0f;
<a name="l01412"></a>01412 
<a name="l01413"></a>01413       <span class="keywordtype">int</span> ishell = my_noscco - 1;     <span class="comment">// First shell will have number 0</span>
<a name="l01414"></a>01414       <span class="keywordflow">for</span> (i__ = 0; i__ &lt; (my_noscco-1); i__++)    <span class="comment">// !!DeBuG!! Iterate to (my_noscco-1) only: the last oscillator is excited in case all other fail (no point in double checking) ??</span>
<a name="l01415"></a>01415       {
<a name="l01416"></a>01416         pac += cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a181131326c64f04da9eb7f199604371e">fco</a>[*mat + i__*<a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * rn[i__];   <span class="comment">// !!DeBuG!! pac[] is calculated on the fly to save registers!</span>
<a name="l01417"></a>01417         <span class="keywordflow">if</span> (pac &gt; temp)       <span class="comment">//  pac[] is calculated on the fly to save registers!  </span>
<a name="l01418"></a>01418         {
<a name="l01419"></a>01419             ishell = i__;
<a name="l01420"></a>01420             <span class="keywordflow">break</span>;
<a name="l01421"></a>01421         }
<a name="l01422"></a>01422       }
<a name="l01423"></a>01423 
<a name="l01424"></a>01424     <span class="comment">/*  ****  Projected momentum of the target electron. */</span>
<a name="l01425"></a>01425       temp = <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed) * rn[ishell];
<a name="l01426"></a>01426 
<a name="l01427"></a>01427       <span class="keywordflow">if</span> (temp &lt; 0.5f)
<a name="l01428"></a>01428       {
<a name="l01429"></a>01429         pzomc = (0.70710678118654502f - sqrtf(0.5f - logf(temp + temp))) / (cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat + ishell * <a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * 1.4142135623731f);
<a name="l01430"></a>01430       }
<a name="l01431"></a>01431       <span class="keywordflow">else</span>
<a name="l01432"></a>01432       {
<a name="l01433"></a>01433         pzomc = (sqrtf(0.5f - logf(2.0f - 2.0f*temp)) - 0.70710678118654502f) / (cgco_SHARED-&gt;<a class="code" href="structcompton__struct.html#a2acdfc304123cea06abea723b5b28d38">fj0</a>[*mat + ishell * <a class="code" href="MC-GPU__v1_83_8h.html#a7aa308231ff3dcfd43ac19023e2e6e15" title="Constants values for the Compton and Rayleigh models:">MAX_MATERIALS</a>] * 1.4142135623731f);
<a name="l01434"></a>01434       }
<a name="l01435"></a>01435       <span class="keywordflow">if</span> (pzomc &lt; -1.0f)
<a name="l01436"></a>01436       {
<a name="l01437"></a>01437         <span class="keywordflow">continue</span>;      <span class="comment">// re-start the loop</span>
<a name="l01438"></a>01438       }
<a name="l01439"></a>01439 
<a name="l01440"></a>01440   <span class="comment">/*  ****  F(EP) rejection. */</span>
<a name="l01441"></a>01441       temp = tau * (tau - (*costh_Compton) * 2.f) + 1.f;       <span class="comment">// this variable was originally called &quot;xqc&quot;</span>
<a name="l01442"></a>01442       
<a name="l01443"></a>01443         <span class="comment">// af = sqrt( max_value(temp,1.0e-30f) ) * (tau * (tau - *costh_Compton) / max_value(temp,1.0e-30f) + 1.f);  //!!DeBuG!! Make sure the SQRT argument is never &lt;0, and that I don&#39;t divide by zero!!</span>
<a name="l01444"></a>01444 
<a name="l01445"></a>01445       <span class="keywordflow">if</span> (temp&gt;1.0e-20f)   <span class="comment">// !!DeBuG!! Make sure the SQRT argument is never &lt;0, and that I don&#39;t divide by zero!!</span>
<a name="l01446"></a>01446       {
<a name="l01447"></a>01447         af = sqrtf(temp) * (tau * (tau - ((float)(*costh_Compton))) / temp + 1.f);
<a name="l01448"></a>01448       }
<a name="l01449"></a>01449       <span class="keywordflow">else</span>
<a name="l01450"></a>01450       {
<a name="l01451"></a>01451         <span class="comment">// When using single precision, it is possible (but very uncommon) to get costh_Compton==1 and tau==1; then temp is 0 and &#39;af&#39; can not be calculated (0/0 -&gt; nan). Analysing the results obtained using double precision, we found that &#39;af&#39; would be almost 0 in this situation, with an &quot;average&quot; about ~0.002 (this is just a rough estimation, but using af=0 the value would never be rejected below).</span>
<a name="l01452"></a>01452 
<a name="l01453"></a>01453         af = 0.00200f;    <span class="comment">// !!DeBuG!!</span>
<a name="l01454"></a>01454                 
<a name="l01455"></a>01455 <span class="preprocessor">        #ifndef USING_CUDA</span>
<a name="l01456"></a>01456 <span class="preprocessor"></span>        <span class="keywordflow">if</span> (warning_flag_2&lt;0)
<a name="l01457"></a>01457         {
<a name="l01458"></a>01458             warning_flag_2 = +1;  <span class="comment">// Disable warning, do not show again</span>
<a name="l01459"></a>01459             printf(<span class="stringliteral">&quot;          [... Small numerical precision error detected computing \&quot;af\&quot; in GCOa (this warning will not be repeated)].\n               xqc=%.14f, af(forced)=%.14f, tau=%.14f, costh_Compton=%.14lf\n&quot;</span>, temp, af, tau, *costh_Compton);    <span class="comment">// !!DeBuG!!</span>
<a name="l01460"></a>01460         }
<a name="l01461"></a>01461 <span class="preprocessor">        #endif</span>
<a name="l01462"></a>01462 <span class="preprocessor"></span>      }
<a name="l01463"></a>01463 
<a name="l01464"></a>01464       <span class="keywordflow">if</span> (af &gt; 0.0f)
<a name="l01465"></a>01465       {
<a name="l01466"></a>01466         temp = af * 0.2f + 1.f;    <span class="comment">// this variable was originally called &quot;fpzmax&quot;</span>
<a name="l01467"></a>01467       }
<a name="l01468"></a>01468       <span class="keywordflow">else</span>
<a name="l01469"></a>01469       {
<a name="l01470"></a>01470         temp = 1.f - af * 0.2f;
<a name="l01471"></a>01471       }
<a name="l01472"></a>01472       
<a name="l01473"></a>01473       <span class="keywordflow">if</span> ( <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a3b57f2faa5765ceba5d352641ca72a9f" title="Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version)...">ranecu</a>(seed)*temp &lt; <span class="comment">/*fpz =*/</span>(af * <a class="code" href="MC-GPU__v1_83_8h.html#a1b7b63b09e6ac252980467701bd80cb7" title="Preprocessor macro to calculate maximum and minimum values:">max_value</a>( <a class="code" href="MC-GPU__v1_83_8h.html#a9fa6cc963655f85859604988fa68a9e2">min_value</a>(pzomc,0.2f) , -0.2f ) + 1.f) )
<a name="l01474"></a>01474       {
<a name="l01475"></a>01475         <span class="keywordflow">break</span>;
<a name="l01476"></a>01476       }
<a name="l01477"></a>01477 
<a name="l01478"></a>01478     }
<a name="l01479"></a>01479 
<a name="l01480"></a>01480 <span class="comment">/*  ****  Energy of the scattered photon. */</span>
<a name="l01481"></a>01481     {
<a name="l01482"></a>01482       <span class="keyword">register</span> <span class="keywordtype">float</span> t, b1, b2, temp;
<a name="l01483"></a>01483       t = pzomc * pzomc;
<a name="l01484"></a>01484       b1 = 1.f - t * tau * tau;
<a name="l01485"></a>01485       b2 = 1.f - t * tau * ((float)(*costh_Compton));
<a name="l01486"></a>01486 
<a name="l01487"></a>01487       temp = sqrtf( fabsf(b2 * b2 - b1 * (1.0f - t)) );
<a name="l01488"></a>01488       
<a name="l01489"></a>01489           
<a name="l01490"></a>01490       <span class="keywordflow">if</span> (pzomc &lt; 0.0f)
<a name="l01491"></a>01491          temp *= -1.0f;
<a name="l01492"></a>01492 
<a name="l01493"></a>01493       <span class="comment">// !Error! energy may increase (slightly) due to inacurate calculation!  !!DeBuG!!</span>
<a name="l01494"></a>01494       t = (tau / b1) * (b2 + temp);
<a name="l01495"></a>01495       <span class="keywordflow">if</span> (t &gt; 1.0f)
<a name="l01496"></a>01496       {
<a name="l01497"></a>01497 <span class="preprocessor">        #ifndef USING_CUDA</span>
<a name="l01498"></a>01498 <span class="preprocessor"></span>
<a name="l01499"></a>01499 <span class="preprocessor">        #endif      </span>
<a name="l01500"></a>01500 <span class="preprocessor"></span><span class="preprocessor">        #ifndef USING_CUDA</span>
<a name="l01501"></a>01501 <span class="preprocessor"></span>        <span class="keywordflow">if</span> (warning_flag_3&lt;0)
<a name="l01502"></a>01502         {
<a name="l01503"></a>01503             warning_flag_3 = +1;  <span class="comment">// Disable warning, do not show again</span>
<a name="l01504"></a>01504             printf(<span class="stringliteral">&quot;\n          [... a Compton event tried to increase the x ray energy due to precision error. Keeping initial energy. (This warning will not be repeated.)]\n               scaling=%.14f, costh_Compton=%.14lf\n&quot;</span>, t, *costh_Compton);   <span class="comment">// !!DeBuG!!</span>
<a name="l01505"></a>01505         }
<a name="l01506"></a>01506 <span class="preprocessor">        #endif</span>
<a name="l01507"></a>01507 <span class="preprocessor"></span>        
<a name="l01508"></a>01508         t = 1.0f; <span class="comment">// !!DeBuG!! Avoid increasing energy by hand!!! not nice!!</span>
<a name="l01509"></a>01509       }
<a name="l01510"></a>01510 
<a name="l01511"></a>01511       (*energy) *= t;
<a name="l01512"></a>01512        <span class="comment">// (*energy) *= (tau / b1) * (b2 + temp);    //  Original PENELOPE code</span>
<a name="l01513"></a>01513     }
<a name="l01514"></a>01514     
<a name="l01515"></a>01515 }  <span class="comment">// [End subroutine GCOa]</span>
<a name="l01516"></a>01516 
<a name="l01517"></a>01517 
<a name="l01518"></a>01518 <span class="comment"></span>
<a name="l01519"></a>01519 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l01520"></a>01520 <span class="comment"></span>
<a name="l01521"></a>01521 <span class="comment"></span>
<a name="l01522"></a>01522 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l01523"></a>01523 <span class="comment">//!  Tally the depose deposited inside each material.</span>
<a name="l01524"></a>01524 <span class="comment">//!  This function is called whenever a particle suffers a Compton or photoelectric</span>
<a name="l01525"></a>01525 <span class="comment">//!  interaction. The energy released in each interaction is added and later in the </span>
<a name="l01526"></a>01526 <span class="comment">//!  report function the total deposited energy is divided by the total mass of the </span>
<a name="l01527"></a>01527 <span class="comment">//!  material in the voxelized object to get the dose. This naturally accounts for</span>
<a name="l01528"></a>01528 <span class="comment">//!  multiple densities for voxels with the same material (not all voxels have same mass).</span>
<a name="l01529"></a>01529 <span class="comment">//!  Electrons are not transported in MC-GPU and therefore we are approximating</span>
<a name="l01530"></a>01530 <span class="comment">//!  that the dose is equal to the KERMA (energy released by the photons alone).</span>
<a name="l01531"></a>01531 <span class="comment">//!  This approximation is acceptable when there is electronic equilibrium and </span>
<a name="l01532"></a>01532 <span class="comment">//!  when the range of the secondary electrons is shorter than the organ size. </span>
<a name="l01533"></a>01533 <span class="comment">//!</span>
<a name="l01534"></a>01534 <span class="comment">//!  The function uses atomic functions for a thread-safe access to the GPU memory.</span>
<a name="l01535"></a>01535 <span class="comment">//!  We can check if this tally was disabled in the input file checking if the array</span>
<a name="l01536"></a>01536 <span class="comment">//!  materials_dose was allocated in the GPU (disabled if pointer = NULL).</span>
<a name="l01537"></a>01537 <span class="comment">//!</span>
<a name="l01538"></a>01538 <span class="comment">//!</span>
<a name="l01539"></a>01539 <span class="comment">//!       @param[in] Edep   Energy deposited in the interaction</span>
<a name="l01540"></a>01540 <span class="comment">//!       @param[in] material   Current material id number</span>
<a name="l01541"></a>01541 <span class="comment">//!       @param[out] materials_dose   ulonglong2 array storing the mateials dose [in eV/g] and dose^2 (ie, uncertainty).</span>
<a name="l01542"></a>01542 <span class="comment">////////////////////////////////////////////////////////////////////////////////</span>
<a name="l01543"></a>01543 <span class="comment"></span><span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01544"></a>01544 <span class="preprocessor"></span>__device__
<a name="l01545"></a>01545 <span class="preprocessor">#endif</span>
<a name="l01546"></a>01546 <span class="preprocessor"></span><span class="keyword">inline</span> 
<a name="l01547"></a><a class="code" href="MC-GPU__v1_83_8h.html#a45a0719009ba55cca7af72f08f1090fb">01547</a> <span class="keywordtype">void</span> <a class="code" href="MC-GPU__kernel__v1_83_8cu.html#a45a0719009ba55cca7af72f08f1090fb" title="Tally the depose deposited inside each material.">tally_materials_dose</a>(<span class="keywordtype">float</span>* Edep, <span class="keywordtype">int</span>* material, <a class="code" href="structulonglong2.html">ulonglong2</a>* materials_dose)
<a name="l01548"></a>01548 {
<a name="l01549"></a>01549       
<a name="l01550"></a>01550 <span class="comment">// !!DeBuG!! The energy can be tallied directly with atomicAdd in global memory or using shared memory first and then global for whole block if too slow. With the initial testing it looks like using global memory is already very fast!</span>
<a name="l01551"></a>01551 
<a name="l01552"></a>01552 <span class="comment">// !!DeBuG!! WARNING: with many histories and few materials the materials_dose integer variables may overflow!! Using double precision floats would be better. Single precision is not good enough because adding small energies to a large counter would give problems.</span>
<a name="l01553"></a>01553 
<a name="l01554"></a>01554 <span class="preprocessor">#ifdef USING_CUDA</span>
<a name="l01555"></a>01555 <span class="preprocessor"></span>  atomicAdd(&amp;materials_dose[*material].x, __float2ull_rn((*Edep)*<a class="code" href="MC-GPU__v1_83_8h.html#ac0c03cf5e803404f5dc2a42467395a3a" title="Value to scale the deposited energy in the pixels so that it can be stored as a long long integer ins...">SCALE_eV</a>) );  <span class="comment">// Energy deposited at the material, scaled by the factor SCALE_eV and rounded.</span>
<a name="l01556"></a>01556   atomicAdd(&amp;materials_dose[*material].y, __float2ull_rn((*Edep)*(*Edep)) );   <span class="comment">// Square of the dose to estimate standard deviation (not using SCALE_eV for std_dev to prevent overflow)</span>
<a name="l01557"></a>01557 <span class="preprocessor">#else</span>
<a name="l01558"></a>01558 <span class="preprocessor"></span>  materials_dose[*material].<a class="code" href="structulonglong2.html#a0d94e96b26765bf28a1a4b76391b88bb">x</a> += (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)((*Edep)*SCALE_eV + 0.5f);
<a name="l01559"></a>01559   materials_dose[*material].<a class="code" href="structulonglong2.html#a8453eaa1efda22c449296aa2ebcf0c7f">y</a> += (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span> <span class="keywordtype">long</span> int)((*Edep)*(*Edep) + 0.5f);
<a name="l01560"></a>01560 <span class="preprocessor">#endif     </span>
<a name="l01561"></a>01561 <span class="preprocessor"></span>          
<a name="l01562"></a>01562   <span class="keywordflow">return</span>;
<a name="l01563"></a>01563 }
</pre></div></div><!-- contents -->
</div>
  <div id="nav-path" class="navpath">
    <ul>
      <li class="navelem"><a class="el" href="MC-GPU__kernel__v1_83_8cu.html">MC-GPU_kernel_v1.3.cu</a>      </li>

    <li class="footer">Generated on Mon Dec 17 2012 13:26:39 for MC-GPU by
    <a href="http://www.doxygen.org/index.html">
    <img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.6.1 </li>
   </ul>
 </div>


</body>
</html>
